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PREFACE 


This brochure contains some 30 programs written for a 
specific programmable pocket calculator, the HP-25. 
These programs implement algorithms in number theory, 
equation solving, algebraic stability theory, numeri- 
cal integration, as well as for the evaluation of 
special functions such as the gamma function, various 
Bessel functions, and the Riemann zeta function. By 
means of the flow diagrams and the detailed descrip-— 
tions that are provided, the programs are easily 
adapted to run on any calculator of comparable capa- 
city. 

The purpose to be served by these programs is di- 
dactic as well as practical. 

The immediate didactic purpose is to enable stu- 
dents of numerical and computational analysis to gain 
first-hand experience with some modern techniques of 
scientific computing. I have used a number of the 
programs in the classroom to demonstrate the numeri-~ 
cal performance of such techniques, using data provi- 
ded spontaneously by the students. Such demonstra- 
tions greatly increase the practical, “know-how" con- 
tent of numerical analysis courses. They also elimi- 
nate the suspicion that the instructor uses examples 
that are rigged to show an algorithm in an especially 
favorable light. 

My purpose is didactic in another, potentially 
even more important sense. Programmable calculators 
have now been in use for almost a third of a century. 
Yet many members of the scientific community, al- 
though they use mathematics in their work, have re- 
mained ignorant of the essentials of modern scienti- 
fic computation. This even holds for some mathemati- 
cians, if they did not choose to specialize in nume- 
rical analysis or computer science. All this is 
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destined to change with the advent of the program- 
mable pocket calculator. By providing the basic fea- 
tures of branching and iteration, these amazing in- 
struments put automatic computation within the reach, 
in the privacy of one's own study, of anybody who re- 
members his basic calculus._No special programming 
language must be learned to operate oné of these de- 
vices; to itS Owner, the access is immediate, and 
there are no waiting times beyond those required by 
the computation itself. An error in the program or in 
the data_can be corrected immediately.-Thus these - 
calculators are interactive to a degree that up to 
now was enjoyed only by the privileged few who were 
connected to an expensive computer by their own ter- 
minal. By publishing my programs I hope to bring some 
of the flavor of modern automatic computation to 
those who so far have been untouched by it. 

Finally, some of my programs, especially those 
dealing with special functions, also serve an emi- 
nently practical purpose. If I quickly need some va- 
lues of a higher transcendental function such as a 
Bessel function of nonintegral order, I know of no 
way to get them with less fuss than by using my 
little calculator. Similar programs could (and un- 
doubtedly will) be written for functions other than 
those represented in my collection; my programs mere- 
ly provide a small sample of what can be done in this 
area. 

Coming back to the didactic aspect, it must be 
admitted that not all areas of modern numerical com- 
putation can be demonstrated satisfactorily ona 
pocket calculator. Because of the limited memory ca- 
pacity, problems dealing with large sets of data, 
such as computations in linear algebra, optimization, 
or partial differential equations, cannot be solved 
on a pocket calculator as they would be solved ona 
large computer. Even in problems requiring small sets 
of data, the lack of memory space imposes limitations. 
Thus it happens, for instance, that my programs for 
determining the zeros or the stability properties of 
polynomials in most cases can deal only with polyno- 
mials of degree not exceeding four. Such programs ne- 


vertheless can serve as models - or "pilot programs" 
- for similar programs designéd to deal with larger 
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sets of data on correspondingly larger computers. In- 
deed, some of my programs for the evaluation of spe- 
cial functions would be written in much the same way 
also for the very largest computers; the only. diffa-—- 
rence here lies in the speed of execution. 

Naturally, this collection is not meant to be ex- 
haustive. The creative student can find ample room 
for the design of new programs, particularly in the 
areas of equation solving, numerical integration, and 
the evaluation of special functions. Furthermore, it 
is likely that many of the programs presented here 
could be improved in some way, say, by making them 
faster, easier to use, or more automatic. I will be 
genuinely pleased to learn about such improvements. 

Already when compiling these programs I have be- 
nefitted from the advice and help of numerous colle- 
gues, students, and friends. J. Waldvogel and D. D. 
Warner contributed ideas. W. Seewald and A. St&hli 
suggested improvements in existing programs. 

E. Specker offered a solution to the problem of wri- 
ting viable power series programs. My wife, Marie- 
Louise Henrici, did most of the checking and proof- 
reading and eliminated numerous errors. The reader 
may judge for himself the excellent work of Brigitte 
Knecht who transformed my manuscript into camera- 
ready copy. The staff of John Wiley & Sons disposed 
of all problems of editing and production with pro- 
fessional know-how. To all these individuals, and to 
any other helpers who may not have been named here, 
I offer my heartfelt thanks. 


Peter Henrici 
Zurich, Switzerland 
February 1977 
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INTRODUCTION 
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This is not a textbook on the programming of pocket 
calculators. To learn how to operate your calculator, 
consult your manual. As a first introduction to the 
manifold uses of the pocket calculator, the recent 
book by Jon M. Smith, "Scientific analysis on the 
pocket ealculator® (Witey, 1975) cai ba cecommendeds 
7 Each program description is arranged according 


to the following scheme: 


1. Purpose, 

2. Method, 

3. Flow diagram, 

4. Storage and program, 

5. Operating instructions, 
6 


. Examples and timing. 


The statement of purpose usually is very brief. 

It is intended to state, in language as nontechnical 
as possible, what the program is good for. 

The description of the method is strictly confined 
to what is being done. For the analytical justifica- 
tion and for technical details, the reader is referred 
to the literature. References to books by the author 


are given in the following abbreviated form: 


4 Introduction 


ENA = "Elements of numerical analysis", 
Wiley, New York, 1964, 
ACCA I = "Applied and computational complex ana- 


lysis", Vol. I, Wiley, New York, 1974, 
ACCA II = "Applied and computational complex ana~ 
lysis", Vol. II, Wiley, New York, 1977. 


The flow diagrams do not follow any particular 
standard format. The intent is simply to make the 
structure of a program visible at a glance. In some 
cases explanations are given for the technical de- 
tails of a program or for abbreviations that are 
used in the flow diagram. 

The section entitled storage and program, natural- 
ly, is the core part of each program description. We 
use the symbol Ry, to denote the storage register num- 
bered k. Quantities that must be stored in the storage 


registers ahead of the computation are encased, thus: 


The program listing uses the key symbols as they are 
found on the keys of your HP-25 calculator, with the 
exception of multiplication, which is indicated by * 
in order to avoid confusion with the letter "x". The 
listing is as compact as possible. For yellow and 
blue instructions it is understood that the keys "f" 
and "g" are to be pressed in advance; no ambiguity 


can arise by our abbreviated notation. The arrows + 
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in front of certain instruction numbers are not part 
of the program. They simply indicate an entrance from 
a GTO instruction. For those readers who would have 
liked to see more complete descriptions of the tech- 
nical details of a program, it must be said that such 
paraphrases can make for extremely tedious reading. 
The best way to understand the details of a program 
is for the user to try to write his own. He then will 
often realize why the particular limitations of the 
pocket calculator forced the author to deviate from 
the more straightforward ways in which such programs 
would be written for a larger computer. 

An attempt has been made to keep the operating 
instructions self-contained, except in a few obvious 
cases. Here we also give indications and explanations 
of possible failures of some programs. 

The examples are not intended to be particularly 
sophisticated. The user should try one of the simpler 
examples to see whether he has correctly loaded the 
program and followed the instructions. In some cases 
the examples indicate how the program would be used 
in the classroom, illustrating particular points of 


analysis. All timings given are approximate. 
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PRIME FACTOR DECOMPOSITION 


1. Purpose 


To produce, for any positive integer m < 10°, its de- 


composition in prime factors. Factors are to appear in 


increasing order. 


2. Method 
We divide m successively by the numbers 
a = 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, ... , 


that is, by 2, 3 and by all odd integers > 3 that are 
not divisible by 3, until either (a) d > Ym or (b) the 
remainder is zero. In case (a), m is prime; the calcu- 
lator shows m and then zero to indicate that the de- 
composition is terminated. In case (b), dis a divisor 
of m. The calculator shows d, replaces m by m/d, and 
starts again, using the last d as its first divisor. 
To compute the d, we start with d = 2 and calcu- 


late the successor d' of d by the formula 
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2-1, if a4 
d+3+e¢, if a>4 


Here the value of € alternates between +1 and -1l, be- 


ginning with -l. 
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3. Flow Diagram 
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4. Storage and Program 


Ry Ry R, R, Ry Re Re 
m vm d € 

00 25 1 
01 sto 0 26 STO-2 
02 2 27 GTO 09 
03 STO 2 + 28 3 
04 1 29 sTO+2 
05 STO 3 30 RCL 3 
06 RCL O 31 CHS 

+ 07 vx 32 «STO 3 
08 sTOlL 33. STO+2 

+ 09 RCL 2 34 GTO 09 
10 RCL1 + 35 RCL O 
ll x<y 36 RCL 2 
12 GTO 41 37 R/S 
13. RCL O 38 + 
14 RCL 2 39 sTOo 0 
15 + 40 GTO 07 
16 FRAC + 41 RCL 0 
17 x = 0 42 R/S 
18 GTO 35 43 0 
19 4 44 GTO 00 
20 RCL 2 45 
21 xZy 46 
22 GTO 28 47 
23 2 48 
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5. Operating Instructions 


After loading the program, move the operating switch 


to RUN. Then press 


FIX 0 


to get the nine-digit display of all integers. If 
prime factors of m are desired, load m into the X 


register. When 


PRGM 
R/S 


is pressed, the calculation will start, then stop by 


displaying smallest prime factor of m. When 


R/S 


is pressed, the next prime factor will be displayed, 
and so forth, until 0 is shown, indicating that all 


prime factors have been found. 


6. Examples and Timing 


36 = 2 * 2 * 3 * 3. Time required about 8 sec. 
71489 = 11 * 67 * 97. Time required about 28 sec. 
987654321 = 3 * 3 * 17 * 17 * 379721. Time re- 


quired about 220 sec. 


ERIE 
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EUCLIDEAN ALGORITHM 


1. Purpose 


To determine the greatest common divisor of two inte- 


gers a and b, not both zero, fal, |b| < 107. 


2. Method 


The Euclidean algorithm. We determine a sequence of 


nonnegative integers {nj} by 


3 
"1 


o t= max (jal,|bl)> , 


n, := min (lal,|bl) , 
eer eee Fedt 
k k-2 7 [nyo |Pk-1 


k = 2, 3, ... . This sequence is decreasing, and thus 


it eventually reaches zero. If ny = 0, then 
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is the greatest common divisor of a and b (see ACCA II, 


§ 12.2). 


3. Flow Diagram 


Display 
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4. Storage and Program 
R R R R R R R 


0 1 2 3 4 5 6 
[= 3) 
No ny 

00 25. GTO 00 
01 ABS 26 
02 xy 27 
03 ABS 28 
04 x <y 29 
05 xy 30 
06 STO 0 31 
07 R¢¥ 32 
08 sTOol 33 

+09 RCL 1 34 
10 x=0 35 
11 GTO 24 36 
12 RCL O 37 
13. RCL O 38 
14. RCL 1 39 
15 + 40 
16s INT 41 
17. RCL 1 42 
18 STO 0 43 
19 * 44 
20 = 45 
21 sto l 46 
22 PAUSE 47 
23. GTO 09 48 


Euclidean Algorithm 17 


5. Operating Instructions 


Load the program. Switch to RUN. Press 
FIX 0 


to get integers displayed as integers. Load a into 
R, and b into R,. Press 


0 1 
PRGM 
R/S 


to start computation. The calculator will stop by 


displaying c, the greatest common divisor of a and b. 
If both a and b are zero, c = 0. 


6. Examples and Timing 


a = 45, b = 96 yields c = 3. Computing time 
about 3 sec. 

a = ~ 965302379, b = 980051 yields c = 997, the 
correct result. Computing time about 5 sec. 

[3] a=1, b= 0 yields c = 1. 
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RATIONAL BINOMIAL COEFFICIENTS 


1. Purpose 


, 


To represent, for arbitrary rational p = E the bino- 


mial coefficients 


z= (8) ea 
bole) = Cy) tel, 


we Py .- Plo-l) (p-2)... (p-ntl) 
BA) 


n=l1, 2, ..., aS irreducible rational fractions, 
Th 
bite) Fa (1) 
n 
where ry and s. are integers, (rs) = 1, and 


Ss) > 0. Such rational values of the binomial coeffi- 


cients are often needed in analytical work involving 
the binomial series, 


(l+x)? = b,(p) x shee. <7 
n=0 


Rational Bhromlad Cae ptieie se: 


2. Method 


The recurrence relation by (nr) \, 


= oon - . 
Bog (9) ntl bf) 7 in O, ly 4y wae , 
is used. If p is a rational number, j) = ae Be ae Be 
is represented in the form (1), this yielan 


xr cr 
n+l P- qn on 


s 
n+l 


q (n+l) sy 


We thus first compute the integers 


r* := (p - qn)r. ’ She+l := q(ntljs_ . 


then by means of the Euclidean algorithm determize 


their greatest common divisor 


c := (r* x 
( n+l’ “n+l 


and find 
x * 
2 - Fn+l 4 _ ntl 
n+1 c ‘ n+1 [ec] . 
The process is started with r, = s, = 1. The short 


0 0 
version of the Euclidean algorithm used here is bor- 


rowed from the program "Exact continued fractions for 


20 Number Theory 


quadratic irrationalities". 


3. Flow Diagram 


(r*,s*) 


(Euclidean algorithm) 
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4. Storage and Program 


Ry Ry R, R, Ry Ky. ve VM, 
00 25  $sTo 7 
Ol CLX > 26 RCL 6 
02 STO 4 27 RCL 6 
03 1 28 RCL 7 
04 STO 2 29 > 
0S STO 3 30 INT 

> 06 RCL 2 31 RCL 
07 PAUSE 32 STO 
08 RCL 3 33 i 
09 R/S 34 7 
10 RCL 0 35 STO 
11 RCL 1 36 x #0 
12 RCL 4 37 GTO 26 
13 * 38 RCL 6 
14 - 39 ABS 
15 STO*2 40 STO=2 
16 1 41 STOs3 
17 STO+4 42 GTO 06 
18 RCL 1 43 
19 RCL 4 44 
20 * 45 
21 STO*3 46 
22 RCL 2 47 
23 STO 6 48 


24 RCL 3 49 
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5. Operating Instructions 


Load the program; move operating switch to RUN. Press 


FIX 
0 


to obtain integer display. Load numerator p of p into 


R, and denominator g into R On pressing 


L° 
PRGM 
R/S 


calculator briefly displays ty) = l and halts by dis- 


playing Sy = 1. Pressing 
R/S 


will cause display of r, and, after a short pause, 


1 
Sy" etc. Pressing 


after display of Sy will cause the decimal value of 


vn 
bn(p) = = 


=) 


to be displayed, without disturbing subsequent compu- 
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tations. Also, after display of Sue the current value 


of n may be exhibited by pressing 
RCL 4 


To re-inspect oe after display of sy! press 


> 
< y 


For large n, integer overflow may cause the values of 


rh and sy to be inaccurate; however, the decimal va- 


lues of bo (pe) will still be accurate. 


6. Examples and timing 


For p = 6, i.e. p = 6, g = 1, we get 
@) = 1, 6, 15, 20, 15, 6, 1, 0, 0, 0, ‘ 
p= 54 ine p = 1, q = 2 yields the values 
UA abs be Sede ely eee te pel 
n n2* 8' 16” 128’ 256’ 1024’ 
33 _ 429 715 _ 2431 4199 
2048’ 32768’ 65536’ 262144’ 524288 
_ 29393 52003 — _185725 
4194304’ 8388608’ 33554432’ 
334305 9694845 20036013 


67108864’ 2147483648’ 4867629602° 


Computing time approx. 165 sec. 
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CONTINUED FRACTION REPRESENTATIONS OF REAL NUMBERS 


1. Purpose 


Given a real number p, to find its standard represen- 


tation either by a terminating continued fraction, 


_ 1 1 1 
p= bp + 5 rite. b 3 (1) 
1 2 n 


if p is rational, or by a nonterminating continued 


fraction, 


+ Fe caver ae ag (2) 


if p is irrational. Here Bo is an integer, and the 


b, (i =1, 2, ... ), if defined, are positive integers. 


2. Method 


By denoting by [x] the greatest integer not exceeding 
x, the b, may be computed as follows (see ACCA II, 
§ 12.2): Let 
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by t= [e] , Pp, t= pa~ by, 
and for k = 1, 2, ... , if Py x 0, 
1 1 

b, := [=] r Pp 1S = Be (3) 

k Py k+1 Py k 
If p is an integer, then already P, = 0, and the re- 
presentation (1) reduces to po = by- If p is rational 
but not an integer, then eventually Oa = 0 for some 
n> 0. The bor bye saves Sy ba thus generated are the 


partial denominators of the continued fraction (1). 


If p is irrational, then Py %# 0 for all k 2 1. 


3. Flow Diagram 
A slight complication arises through the fact that 
the operation INTp agrees with [p] only if p 2 0, or 


if p < 0 and p is an integer. For nonintegral p < 0 


we have 
[ep] = INTp - 1. 


Otherwise, the flow diagram is straightforward. 
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Display 
b, := INTp - l 


Display 
b, := INT p 


Display 0 


Continued Fraction 


4, Storage and Program 


R 


0 


RF 


s. 


Reo eee it OG bettas 


oA 


s/ 
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5. Operating Instructions 


Load the program; switch to RUN. Load p into the X 


register. Press 
FIX 0 
to get the integer display. By pressing 


PRGM 
R/S 


bo will be displayed. Pressing 


R/S 


repeatedly will display bi, b If a zero is 


be Git ets 
displayed, this means that fe algorithm has termina-~- 
ted, the number p being rational. Note, however, that 
because of rounding errors the algorithm does not al- 
waysS operate as predicted by the theory. Thus for ra- 
tional p we may sometimes get an exceedingly large 

Dit in place of 0, indicating that on the computer 


Prt 
irrational p, the b, from a certain point onward ge- 


was not exactly zero as it should have been. For 


nerally are wrong. These deficiencies could only be 
corrected by working with multiple precision (compli- 
cated on the HP 25) or, in certain special cases, 


with integer arithmetic (see the following program). 


6. 
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Examples 


2/3 yields b, = 0, 1, 2, 0. Indeed, 


D 
Wt 


ze) 
It 
( 
ae 
wn 
Ke 
BE 
0) 
rR 
QR 
1) 
oa 
Il 
i] 
wn 
. 
N 
H 
=] 
Q 
oO 
Qa 


re) 26/47 yields bi = 0, 1, 1, 4, 4, 1, 
(20000000). The last entry is wrong. In fact, 


26 1 i 1 1 i 
4 -ppa le ale eal. 


For p = e we get, as discovered by Euler, b; = 


SI 


The last entry should be 10 and hence is wrong. 
p = 7 furnishes bi = 
3, 7, 15, 1, (293) . 


The first few approximants of the continued frac- 


tion for n thus are 
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+ 3.14285714 , 
1 1 = 333 

vee = ipg = 3-14150943 , 
1 it ie 355 
YX +a + x arg = 3-14159292 . 


The next approximant, unfortunately, is already 


wrong (b, should be 292), which is indicative of 
the great accuracy that is already achieved by 


the foregoing loworder approximants. 


oom 1+ 75 
6 _—_—_ = = = = a 
{6} p 2 . We get by by aes boy 1, an 
then bo = 2, which is wrong. 


p = v7. The sequence of bis is as follows: 
2, 1, 1, 1, 4, l, 1, 1, 4, l, li 1, 4, Ly Ly l, (3) 


The last two examples illustrate the fact that 
the continued fractions representing quadratic 
irrationalities are ultimately periodic. The 


following program shows how to construct them 
without rounding error. 
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EXACT CONTINUED FRACTIONS FOR 
QUADRATIC IRRATIONALITIES 


1. Purpose 


To construct, without rounding error, the standard 
continued fraction representation of a quadratic ir- 


rationality, that is, of a real number 


> = Rtas 


aed (1) 
wherr p, q, m, xr are integers and r is not a square, 

m # 0. The set of quadratic irrationalities coincides 
with the set of irrational solutions of quadratic 
equations ax? + bx + c = 0 with integer coefficients 
a, b, c. By a celebrated theorem of Lagrange (see 
ACCA II, § 12.2) the standard continued fraction re- 
presentation of a quadratic irrationality is ultimate- 
ly periodic. It thus suffices to determine, in addi- 
tion to the nonperiodic part, one period of the con- 


tinued fraction. 
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2. Method 


The same division algorithm as in the preceding algo- 
rithm, but performed in rational arithmetic. If r is 
fixed, the numbers of the form (1) create a field F; 
that is, sums, differences, products, and quotients 
of such numbers can again be represented in the same 
way. This is clear for sums, differences, and pro- 


ducts; for quotients, we may use the identity 


m - mp + mqvr 


voir 
1 


p + qvr q’r - p* 


cach number 9 € F thus is represented by a triple of 
integers [p, q, m]. Many triples represent the same p, 
but among all triples representing p there is a unique 
representation, which we call the reduced representa- 
tion, that is characterized by the fact that p, q, m 
have no common factors and that m > 0. We may use the 


symbol 
ip, a, m) v 
for any representation of p, and 
[P, a, m] Yop 


for the reduced representation. For any extended nu- - 


merical work in F one should use reduced representa- 
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tions, both to prevent overflow and to recognize re- 
peated elements. To this end, the integers p, q, m 
should be divided after each operation in F by their 
greatest common divisor (g.c.d.), say, dad. This is de- 
termined by first calculating the g.c.d. (p,m) of p 
and m by the Euclidean algorithm, and then construc- 
ting d = ((p,m),q) by another application of the Eu- 
clidean algorithm. 


For our present purposes the division algorithm 


is formulated as follows: Let a9 := 9, and for k = 0, 
1, 2, ... compute 
by := [a,] ’ 
a := : 
k4+1 Om by 


There is no possibility of vanishing denominators, be- 
cause the continued fraction is a priori known to be 
infinite. The period of the continued fraction is com- 
=a, = 0 and n > 0, that is, 


ktn k 
as a reduced triple representing a 


pleted as soon as a 


k is repeated. 
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3. Flow Diagram 


Compute 


n~ (p,q, m] 


Display 


b = [a] 


Compute 


= av , ™m 
b- [a] [P. ‘ 4 


Inspect 


Pr, Gy, m 


Display 
(p,m) 


Display 


((p,m) ,q) 


To save programming space, the program computes bo 
correctly only if p > 0. For similar reasons, a 
branching after the second application of the Eucli- 


dean algorithm must be performed manually. 
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4. Storage and Program 


Ro Ry Ry R, R, Ro Re Ro 
xr temp temp temp 

00 25 R+ 

o1 RCL 3 26 STO 2 
02 vx 27. sto 7 
03 RCL 1 28 RCL 0 
04 - 29 STO 6 
05 RCL 0 > 30 RCL 6 
06 + 31 RCL 6 
07 RCL 2 32 RCL 7 
08 7 33 = 

09 INT 34 INT 

10 R/S 35 RCL 7 
ll RCL 2 36 STO 

12 * 37 id 

13 STO-0 38 - 

14 RCL l 39 STO 7 
15 “re 40 x # 0 
16 RCL 3 41 GTO 30 
17 * 42 RCL 6 
18 RCL 0 43 R/S 

19 x? 44 RCL 1 
20 = 45 STO 7 
21 RCL 2 46 GTO 30 
22 STO*1] 47 STO=0 
23 CHS 48 STO+1 
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A short version of the Euclidean algorithm is con- 


tained in instructions 30 through 43. 


5. Operating Instructions 


Load the program, then move the operating switch to 
RUN. Press 


FIX 0 


to get the integer display. Load 


p into Rp 


q into 


R 
m into R 
r into R 


1 
2 
3 


The following instructions should now be carried out 
cyclically. *Press 


PRGM 
R/S 


The computer very shortly will display bo. Pressing 
R/S 


will cause program to compute nonreduced representa-— 
tion [p, q, m] of 1/ (a - by) and to compute and dis- 


Quadratic Irrationalities 


play (p,m). Pressing 
R/S 
once more causes display of dad = ((p,m),q). Now press 


GTO 47 
R/S 


This will compute the reduced representation [p,q,m] 
%y a): The calculator will stop by displaying d. To 
inspect reduced representation, press 


RCL 0 
RCL l 
RCL 2 


This will in turn show p, q, m. Now go back to * to 


compute b, and so forth. The period of the fraction 


1 
is complete if a reduced triple p,q,m is repeated; 
the period begins where the repeated triple occurs 


the first time. 


6. Examples and Timing 


ep = 24 . We obtain 
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b. 1 5 2 3 2 
1 

3 3 

: ~-1 l 1 1 1 
1 

m 17 2 3 2 3 


The repeated triple is encased. Indicating the period 


by a bar, the sequence of the bi thus is 


[2] p= 4 5 v5 Here the results are 

DQ. 0 l 8 6 7 1 1 
i 

Py 4 28 17 15 15 13 

q, 1 -7 7 7 7 7 7 
m. 7 11 4 5 4 19 11 
b, 1 30 1 1 1 7 6 
Py 5 15 15 5 6 13 is 
q; 7 7 7 7 7 7 
m, 20 1 20 11 19 4 5 


Total computing time about 10 min. 


Part 2 
ITERATION 


Iteration 41 


ITERATION 
1. Purpose 


To determine a fixed point of a given function f, 


that is, a solution of the equation 
x = f(x). (1) 


This may also be used to construct solutions of equa- 


tions of the form 
g(x) = 0 (2) 
by seeking fixed points of the function 
£(x) := x +c g(x), 


where c is a suitably chosen constant or function. 


2. Method 


We construct an iteration sequence {x} by choosing 


xX, arbitrarily and forming 


0 


42 Iteration 


fees ce (3) 


If £ maps a closed interval I into itself, if f is 
contracting on I (i.e., 


that |£(x) - fly) | 


if there exists k < 1 such 

k|x - y| for arbitrary x, y € I), 
it can be shown that the equation (1) has precisely 
one solution s in I, and any iteration sequence star- 


ted with an X) € I converges to s (see ENA, § 4.2). 


3. Flow Diagram 


n+1 
x 3= £(x) 


Display x 


The index n merely counts the iteration steps; it is 


not actually needed in the computation. 


4. Storage and Program 


00 
O1 
02 
03 
04 
05 
06 
07 


STO 0 
CLX 
STO 7 
1 
STO+7 
GTO 10 
STO 0 


10 
ll 
12 
13 
14 
15 


Iteration 


PAUSE (R/S) 
GTO 04 


The program to compute f should be in the locations 


10 through 49. The last instruction of this program 
should be GTO 07. 


5. Operating Instructions 


Load the program; 


mode of displaying numbers, 


Load the starting value x 


sing 


FIX 8 


0 


turn the switch to RUN. Select the 


for example, by pressing 


43 


into the X register. Pres- 
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PRGM 
R/S 


will start the computation. The computer will pause 
briefly after computation of each iterate and display 
it. If a stop is desired, instruction 08 should be 
R/S; in this case, 


R/S 


must be pressed after each display. The iteration 


index n is displayed by pressing 
RCL 7 
after display of x. No convergence test is provided 


by this simple program. 


6. Examples and Timing 


f(x) := /2 + x . The program to compute f is 


10 RCL 0 
ll 2 
12 + 
13 Vx 


14 GTO 07 


Starting with 


*o 


13 
*14 
*15 


Iteration 


0, we get 


-41421356 
-84775907 
-96157056 
-99036945 
-99759091 


en 


= 1.99999996 
= 1.99999999 
= 2.00000000 


which, to the number of digits shown, equals s 


= 2. Computing time about 20 sec. 


This equals s 


Again starting with XQ i= 0 we 
x) = 1.00000000 
0.50000000 
0 .66666667 
0 .60000000 
X19 = 0.61803400 
Xoq = 0.61803398 
X51 = 0.61803399 
X59 = 0 .61803399 


t 
NIE 


(f5 - 1), the unique positive 
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solution of 


Computing time about 30 sec. 


f(x) := (x - 1)?. With x, := 3 we get a rapidly 


divergent sequence; with XQ t= 2 we find 
1.0000 

0.0000 

1.0000 

0.0000 , 


x * K x 
m= Wn 
u 


an oscillating sequence. 


[4] £(x) :=e oa X, ?= 0 produces a sequence that, 
although theoretically convergent, does not in 
fact converge. After n = 833 iterations, the 
sequence cycles between 


0.99954403 and 
0.00004561 


*933 
*934 
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ITERATION WITH AITKEN ACCELERATION 


1. Purpose 

To determine a fixed point of a given function f more 
rapidly than with ordinary iteration. 

2. Method 

Aitken acceleration. Along with the sequence {x}, 


which is generated as in the preceding program, we 


generate the sequence of accelerated values {x'}, 


where 
(ax)? 
, n 
Xi oFe xX 7) 
a x 
n 
(4 = forward difference operator). If {x} converges, 


then under certain conditions (see ENA, § 4.4) {x1} 


converges to the fixed point more rapidly than {x}. 
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3. 


Iteration 


Flow Diagram 


Display X) I= ¥ 


Compute and display XG 


x := 


4. 


Storage and Program 
Ry Ry R, 
ge *o 

00 

ol CLX 

02 STO 7 

03 RCL 0 

04 STO 5 

05 GTO 39 
+ 06 RCL 7 

07 x # 

08 GTO 15 

09 RCL 

10 STO 1 

ll STO 

12 1 

13 STO+7 

14 GTO 39 
> 15 RCL 6 

16 STO 2 

17 PAUSE 

18 RCL 1 

19 ad 

20 RCL 1 

21 RCL 0 

22 - 

23 STO 3 


Iteration with Aitken Acceleration 


> 


ww rP NY oO Ff 


STO 6 


GTO 


49 
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Program locations 39 through 49 are reserved for the 
program computing f(x). This program should assume x 


in R, and should put f(x) into R,. Only the stack may 


6 
be used for temporary storage. The last instruction 
should be GTO 06. The program for computing the f of 


example is shown above. 


5. Operating Instructions 


Load the program (including the program to compute f); 
switch to RUN. Select the mode of displaying numbers, 


for instance 
FIX 8 


If f contains trigonometric functions with argument 
in radians, press 


Load the starting value Xo into Ro: When 
PRGM 
R/S 


is pressed, the computation starts. The calculator 


will briefly display each X42 and stop at display 
of x: Press 


Iteration with Aitken Acceleration 


R/S 


51 


to start a new cycle. No convergence test is provided 


nor are iterations counted. 


6. 


G) 


Examples and Timing 


E(x) := /2 + x . Sequence {x} , if x) = 0: 
-03942606 
-00208254 
-00012537 
-00000776 
-00000048 
-00000003 
-00000000 


MR OR RRS 
De M- B-wW-N-F-O- 
HT] 
NNNNNN WN 


The fixed point s 2 has now been reached in 
six iterations (which requires the computation 
of 8 iterates x,)+ Computing time is about the 
same as before, because of the additional work 


required to compute wm 


f(x) :=e * . with X) = 0 we get 
x, = 0.36787944 x, = 0.61269984 
0.69220063 0.58222610 


0.50047350 0.57170577 
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*16 
A54 


At this 
whereas 


correct 


= 0.60624354 x3 = 0.56863881 
= 0.56706790 X14 = 0.56714330 
= 0.56718605 XI5 = 0.56714329 


point the sequence {x1} has converged, 
the sequence {x} still has only four 


digits. Computing time about 80 sec. 
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AITKEN-STEFFENSEN ITERATION 


1. Purpose 


To compute the fixed points of smooth functions f 


without restrictions on the slope of f. 


2. Method 


Aitken-Steffensen iteration. This is a variant of the 
Aitken acceleration considered in the preceding pro- 


gram. Each Xo is used as a starting point for a new 
iteration. Equivalently, we generate a sequence 


(0) 


{x'™) } by choosing x arbitrarily and forming 


(E(x! ) _ inh? 
(n))) 


x (ntl) = xm) 


(n)) x 


f£ (f(x - 2£(x + 


If f has a fixed point s, if f'(s) # 1, and if (9) 


is chosen sufficiently close to s, then 4 (P) 
quadratic convergence (see ENA, § 4.11). If £'(s) = 1, 
convergence is dubious; in particular, some x i) may 


fail to be defined because of vanishing denominator. 
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3. Flow Diagram 


Display 


-- 


A 


wv 


Cone) 
STS. s 
3,5 


“a 
ty 
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The program to compute f should be in locations 35 
through 49. The value of x is in Ro y = f(x) should 


be put in Re: Use the stack only for temporary storage. 


The last instruction is GTO 06. The program for exam- 
ple 4] is shown above. 


5. Operating Instructions 


Similar to preceding program. After 


PRGM 
R/S 


is pressed, the computation starts. The computer dis- 


plays each x) and stops. Press 


R/S 


to start new cycle. No convergence test is provided. 


6. Examples and Timing 


f(x) := /2 +x. x 69) = 0 yields sequence {x™)} 


as follows: 


es 


= 2.03942606 
x (2) = 2.00000802 
(3) = 2.00000000 
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The limit has been reached (to within machine ac- 


curacy) in three iterations. Computing time 10 


sec. 
= 1 (0) _ ? 
f(x) 25 as = 0 yields 
x?) = 0.66666667 
x'?) = 9.61818182 
x) = 0.61803399 
Computing time 15 sec. 
£(x) s= (x - xy? é gio) = 3 yields 
xf) = 2.75000000 
2.63888889 
2.61864329 
2.61803453 
2 .61803399 


The last value is 5(3 + YS), the larger fixed 


ROI 0 we get 


point of f. Starting with x 
-50000000 
- 38888889 
-38199234 
- 38196601 
-38196601 


oo oO 0 89 


The last value is $(3 - 75), the smaller fixed 
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point. Computing time in both cases about 15 sec. 


[4] £(x) := a te With x9) = 0 we get 


x) 


-50001135 
-32882657 
- 23868709 
-19120788 
~17597093 
- 17456388 
»17455280 
-17455280 


oooooeoe°oeo:9®o 


Iteration converges after eight steps. Computing 
time about 30 sec. 
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NEWTON ITERATION FOR COMPLEX ROOTS 


1. Purpose 


To demonstrate the convergence of Newton's method to 
compute the square root of a given complex number 


c¢ = a+ ib # 0. This method consists in forming the 


sequence {z 0} = (x) + iy} by choosing 2 and then 
computing recursively 
sil C : 
ay 3 (2, + 2 po Oy D2 ad (1) 


This sequence converges to the value of Yc nearest to 
Z if there is precisely one such value. It diverges 


if z, lies on the straight line through O perpendicu- 


0 
lar to the straightline segment joining the two values 


of Ye (see ACCA I, § 6.12). 


2. Method 


We evaluate the sequence (1) until [z - 2 < €, a 


n+1 a 
prescribed tolerance. The number of iteration steps is 


counted. 
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3. Flow Diagram 


Display zy 
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4. Storage and Program 


BG Ry Be Ry By Re Re SF 
a b x y € |z [? Ax n 
0 0 n n 
00 25 RCL 2 
01 CLX 26 * 
02 sto 7 27. - RCL 
+ 03 1 28 RCL 3 
04 STO+7 29 * 
05 RCL 3 30 - 
06 RCL 2 31 RCL 5 
07 >» Pp 32 = 
08 x? 33. RCL 3 
09 «sTO 5 34 - 
10.» RCL 35 2 
lls RCL 36 + 
12 * 37. STO+3 
13. RCL 38 RCL 6 
14 RCL 3 39 sTo+2 
15 i 40 +P 
16 + 41 $TO 6 
17s RCL 5 42 RCL 2 
18 + 43. PAUSE 
19 - RCL 2 44 RCL 4 
20 - 45 RCL 6 
21 2 46 x=y 
22 7 47 GTO 03 
23 +=sTO 6 48 RCL 3 
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5. Operating Instructions 


Load the program. Turn the operating switch to RUN. 


Select the mode of display numbers, for instance 
SCI 8 
Load the data as follows: 
into 


b into 


x into 


nD Dw A 
wn fF Oo 


Yo into 


Load € into Ry- For instance, for € = io? press 


When 


PRGM 
R/S 


is pressed, the calculator starts iterating, pausing 
briefly after each cycle while displaying xe After 


convergence has been achieved, the calculator will 


Newton Iteration for Complex Roots 


stop and display final xo° To find Y,. Press 
x < y 
To display number of iterations required, press 


RCL 7 


If all iterates are to be recorded, 


should be changed to R/S. The calculator will then 
stop after each cycle, displaying xO: To display Yue 


press 
RCL 3 


To compute next iterate, press 


R/S 


6. Examples and Timing 


c= 3+ 4i, z, = 1, ¢ = 107. The successive 


0 
iterates are 


instruction 43 


2.00000000 + 2.00000000 i 
1.87500000 + 1.12500000 i 
1.99632353 + 0.99387255 i 
1.99999968 + 1.00001144 i 
2.00000000 + 1.00000000 i 
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Convergence has been achieved. Computing time 
35 sec. 


c and e€ as before, 


Zp = -i. Algorithm converges 
to 
- 2.00000000 - 1.00000000 i 
inn = 4 iterations. Computing time 40 sec. 


c = -1 + 0.00001 i, z, = 1, € = 10°? 


0 . Algorithm 


converges to 


0.00000500 + 1.00000000 i 
inn= 


24 iterations. Computing time about 2 min. 


[4] ¢=-l, Z) = 1. No convergence. 


Part 3 
POLYNOMIALS 
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HORNER ALGORITHM 


1. Purpose 


Lh be ie 
bhistic, ages ey 


Given an arbitrary polynomial of degree s 4 with real 


coefficients 
= 4 3 2 
p(x) = agx + a,x + a,x + a,x + a, ; 
and given an arbitrary real number Xgr to determine 


the coefficients De in the representation of p in 


powers of h := x - Xoe 


= 4 3 2 
P(X + h) = boh + bjh + bh + b,h + b, ; 


that is, the Taylor coefficients of p at X- 


2. Method 


The Horner algorithm (see ENA § 3.4; ACCA I § 6.1). 


(m) 


We generate coefficients be as follows. Let 


at 
bp! eet SP ee Ws it es, gol, 
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and for m= 0, l, ..., 4: 

(m) ._ ,(m-1). ,(m) |. (m) (m-1) z e. 
by := bo : bo := XgP oi bo ,» ne=tdl,...,4-m. 
Then 

Be DOOM 0 Gy Oy. Dp wees pd 
m m 


3. Flow Diagram 


Because an address modification is not available, the 
programming here is different from what it would be 
on a larger computer. By a cyclic permutation (here 
called rotation) of the ays the coefficients operated 
on are always found in the same registers. After the 
algorithm is executed, the b, are stored where, pre- 


viously, the a. were stored. 


1 
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Display by 


Stop 
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4. Storage and Program 


3 4 5 6 7 
n.m 
0) 
6 
1 
Enashig S- Oar REE 28 sto 0 
puarke ee 04 INT 29 RCL 2 
05 sfo 7 30 STO 1 
+ 06 RCL 7 31 RCL 3 
07 FRAC 32 sTO 2 
08 9 33 RCL 4 
09 * 34 STO 3 
10 RCL 7 35 RCL 6 
ll + 36 STO 4 
12 3 37 . 
13 x <y 38 1 
14 GTO 25 39 STO +7 
sf 15 RCL O zs 40 . 
16 RCL 5 41 5 
17 * 42 RCL 7 
18 = RCL1 43 FRAC 
19 + 44 x <y 
5 20 STO 1 | eee EL 
21 3 46 1 
22 RCL 7 47 STO +7 
23 x2y 48 GTO 03 
24 GTO 49 > 49 RCL O 


Note: A fractional index n.m is used to save storage. 
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5. Operating Instructions 


Load the program and switch to RUN. Load the coeffi- 


cients an into Ro (n = 0, 1, ..., 4) and Xo into Ro - 
When 

PRGM 

R/S 


is pressed, the calculator computes the bo and stops 


by displaying b At the same time, the b are stored 


0° 
in Rn (m = 0, 1, ..., 4), thus enabling the operator 
to continue the computation immediately with a dif- 
ferent XQ- The original a, are lost. 


6. Examples and Timing 


Consider the polynomial 
p(x) = x - 4x” + 3x? - 2x+5. 


Carrying through the algorithm with Xq = 0.5 
yields the Taylor coefficients at 0.5, 


b_ = 1.0000, -2.0000, -1.5000, -1.5000, 4.3125 
(time required about 29 sec). The last coeffi- 


cient is p(0.5). Repeating the algorithm with 
the foregoing data yields the Taylor coefficients 
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at l, 
be = 1.0000, 0.0000, -3.0000, -4.0000, 3.0000 
Repeating three times with x, = -0.33333333 


i) 
yields the coefficient arrays 


1.0000, -1.3333, -2.3333, ~2.1481, 4.0123 
1.0000, -2.6667, -0.3333, -1.1852, 4.5309 
1.0000, -4.0000, 3.0000, -2.0000, 5.0000 


The last array is identical with the initial 
array of coefficients, demonstrating the "group 


property" of the Horner algorithm. 
(2] For the polynomial p(x) =x, X = l we get 


bh = 1.0000, 4.0000, 6.0000, 4.0000, 1.0000, 


the binomial coefficients Cy Computing time 


29 sec, 
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NEWTON'S METHOD FOR POLYNOMIALS 


1. Purpose 


To determine all real zeros of a real polynomial of 
degree 4, 


p(x) = ax +a (1) 


2. Method 


Newton's method, using Horner's scheme to evaluate p 
and p' and deflating the polynomial after each zero 
has been found. Newton's method requires forming the 


sequence {x} according to 


' (2) 


with an arbitrarily chosen starting value Xq- The se- 
quence {x} will converge to a zero, provided a real 
zero exists and x, is chosen close enough to the zero. 


0 


(Unless p'(0) = 0, the choice x, = 0 will frequently 


0 
work and will produce convergence to the zero of 
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smallest modulus. | For the purpose of this algorithm, 
Horner's scheme may be formulated as follows: Compute 


the array {b.} by 


Then p(x) = b, and p'(x) =c The iteration (2) is 


4 
terminated if 


-_ -6 
i = [p(x )| < 100. 


(This tolerance may be changed to any one-digit power 
of 10, but too stiff a tolerance may cause the se- 
quence {x,) to cycle instead of to converge.) If 

p(z) = 0, the b, are simply the coefficients of the 
deflated polynomial 


Because the coefficients by are available, the de- 
flation thus is achieved by moving the bi into the 


locations of the a; (see operating instructions). 
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3. Flow Diagram 


6 


Ip(x) | < 10° 


? 


Compute p'(x) 


Display x 
Stop 


The program presents a particular challenge because 


there are ten essential quantities a a a 


ori 
Dy x that have to be saved. 


203 


ic 
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4. Storage and Program 


) 1 2 
30 Be 3 
A 
= Ol ENTER 
02 ENTER 
03 ENTER 
04 RCL O 
0S * 25 
06 RCL 1 
07 + 
b 08 sto 5 
09 * 
10 = RCL 2 
11 + 
Dy 12 STO 6 
13 * 
14 RCL 3 
15 + 
b, 16 sto 7 
49 * 
18 RCL 4 
PRO 19 + 
20 ABS 
21 EEX 
22 CHS 
23 6 
24 x2y 
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5. Operating Instructions 


Load the program. Turn the operating switch to RUN. 
Select the mode of displaying numbers, for instance, 


by pressing 
FIX 8 


Load the coefficients a, of given polynomial into 
Ris i=0, ..., 4. {Be sure to use indexing of co- 
efficients as in (1) -] Load the starting value Xo 


into the X register. Press 


PRGM 
R/S 


to start computation. The iterates x will be compu- 
ted; each x will be displayed briefly. The calcula- 
tor will stop if x, meets convergence test Ip (x) | < 
Lor displaying final value of xo" Exponent -6 in 
convergence test may be changed to -m by changing in-~ 
struction 23 to m, where m is any one-digit integer. 
To deflate polynomial when zero has been found, 


press 


RCL 0 
sto l 
RCL 5 
STO 2 
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RCL 6 
STO 3 
RCL 7 
STO 4 
CLX 
STO 0 


Load the starting value x, into the X register (if 


0 


x, = 0, it is already there) and press 


0 


R/S 


to restart computation. The process may be repeated 


until all real zeros have been found. 


6. Examples and Timing 


p(x) := ae - 16x? + 72x" - 96x + 24. 


Starting with Xp = 0 we get the smallest zero 


z, = 0.32254769 


after four iterations. Deflating and again star- 


ting with Xp = 0 yields 


zy = 1.74576110 


after five iterations. Deflating and starting 
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with XQ = 0 after five iterations yields 


23 = 4.53662030 


Deflating once more yields 


9.39507092 


N 
tt 


in one iteration. Computing time about 4 sec 
per iteration. 
p(x) <:= - tox + 35x" - 50x + 24 

= (x - 4) (x - 3) (x - 2)(x - 1). 


The program correctly finds the zeros 


z= 1.00000000 
Zz = 2.00000000 
z= 3.00000000 
2, = 4.00000000 


Total elapsed time (including deflations) 
2.75 min. 
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BERNOULLI'S METHOD FOR SINGLE DOMINANT ZERO 


1. Purpose 


To determine the zero of greatest absolute value of a 


real polynomial of degree 4, 


4 
p(x) = x +a,x +a 


2 
1 2 3 4 ' () 


if there is a single such zero. 


2. Method 


Bernoulli's method (see ENA, § 7.1; ACCA I, § 7.4). 


A sequence {x} is generated by choosing starting 


values Xe Xr Kos Xy “arbitrarily” and letting 
Mae. SAA pea SS otnte: P Stina? By kaeg? 
n= 4, 5, ... . One forms the quotients 
*n 
q, = yO (2) 
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if p has a single dominant zero w and if the starting 
values are not in a certain set of measure zero, then 
the q, are defined, at least for sufficiently large n, 


and 


w = lim qn, 


n--e 


In this program the calculation is stopped if 


< 10”. (3) 


This may cause accidental termination if two consecu- 
tive q's happen to be identical. The program will al- 
so halt if an x becomes accidentally zero, or if 
there is overflow or underflow because Xn exceeds the 
range of the computer. 

Experience has shown that for polynomials with 
small integral coefficients these accidents most often 
are caused by a too special choice of the starting va- 


lues, such as x, = x, = x = 1. The program, 


0 1 2 3 
therefore, generates its own highly irrational star- 


= xX 


ting values. 


3. Flow Diagram 


A difficulty is generated by the fact that all avai- 
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lable storage locations are taken up by the four co- 
efficients of the polynomial and by the four currently 
needed XW° Therefore, there is no orderly transition 
from n to n + 1; rather, each x is overwritten with 
X41] 48 soon as it is no longer needed. Furthermore, 
the quotients qn and Gael are recalculated each time 


they are needed. 


Generate Xgs Xy0 X50 a 


yes 


Display w = qd 
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4. Storage and Program 


i) 1 2 3 4 5 
ei 22 ae a4 "5 me 
O 00 25 
ol 1 26 
X, 02 sTo4 F 27 
03 YX 
x, 04 220.2 ee 29 
05 e 2 30 
rs 06 STO 6 31 
07 log 32 
x, 08  SsTO 7 = 33 
O09 RCL Say 34 
10 RCL 7 x, 35 
ll * ay, 36 
12 RCL 2 37 
13. RCL 6 38 
14. STO 7x: x, 39 
15 * 40 
16 + 41 
17) RCL 1 42 
18 RCL 5 43 
19 STO 6 Xt D 44 
20 * 45 
21 + 46 
22 RCL O 47 
23 RCL 
24 sto s kyte ty 49 


48 lL 


6 7 
xy Xo 
* 

+ 
CHS 


28 Sto 4 Xp 


RCL 5 
PAUSE (or R/S) 43 
ENTER 


* ENTER 


RCL 5 *2 
RCL 6 XxX, > 4a 


x < y 


GTO 09 


R¥ 
R ¥ 


GTO 00 
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5. Operating Instructions 


Load the program, switch to RUN, and choose mode of 


displaying numbers desired, for instance, by pressing 
SCI 8 


to get floating eight-digit display. Another possibi- 
lity is initially to display few digits only, which 
makes it easier to check convergence "by eye." Be sure 
to define coefficients of polynomial as in (1). Load 
the coefficient a, into Ry (k = 1, 2, 3, 4). Press 
PRGM 
R/S 


to start computation. The calculator will display 
briefly each qn° (If an indeterminate stop is desired, 
instruction 31 should be changed to R/S. In this case, 
R/S must be pressed after each display to continue 
computation.) The calculation will terminate if con- 
dition (3) is satisfied; the last q,, will be dis- 
played. 

If starting values other than those provided by 
the program are desired, these should be loaded as 
follows: 

into 
into 
into 


x3 into Rq - 
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In this case, the calculation is started by pressing 


GTO 09 
R/S 
6. Examples and Timing 
j ow 5 Ar 
p(x) := x" - 10x” + 35x° - 50x + 24 


(x - 4) (x - 3) (x - 2) (x - 1). 


The calculator displays 


w = 4.0000000 


after 60 iterations. Computing time approximately 


4 min. 


[2] p(x) = x4 + 3x? - 7x? - 15x + 18 
(x + 3)7 (x - 1) (x - 2). 


Convergence is very slow because the dominant 


zero has multiplicity > 1. Faster convergence is 
achieved (see ENA, § 7.4) by choosing the star- 


ting values 


Xp = 7 a, = 3 

x1 = - (2a, + a)Xo) = 23 

Xo = 7 (3a, + ayXy + a,x) = - 45 

x, 57 (4a, + a3Xp + a,x, + a,X) = 179 


This yields 
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w = -2.9999999 


after 45 iterations. Computing time approximately 


3 min. 
p(x) r= x! - x? - 9x” - 2x - 3 yields 
w = 3.6652758 


Computing time 3.4 min. 
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BERNOULLI'S METHOD FOR COMPLEX 
CONJUGATE DOMINANT ZEROS 


1. Purpose 


To determine the quadratic factor formed by the two 
zeros of largest modulus of a real polynomial of de- 
gree 4, 


p(x) = a + a,x” +a “2 +t+a.xta,, (1) 


provided that all remaining zeros have smaller modulus. 


2. Method 


Bernoulli's method (see ENA, § 7.5), which in this 

situation is identical with a pristine version of the 
quotient-difference algorithm (see ACCA I, § 7.6). We 
generate the sequence {x)} as solution of the diffe- 


rence equation 


xy = - (ax iat a,x + 
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n = 4, using arbitrarily chosen initial values Xqe 


Xor %3- With the x, we form the quotients 


x 


1’ 


as in the preceding algorithm, and with the quotients 


the differences 


Finally, we form the quantities 
n-1 * 4n-1 ' 8p ** Ip-29n-1 * 


If p has two dominant zeros Wy and Woe and if all 
other zeros have smaller modulus, then the limits 


’ s := lims 
n--o noo 


exist and satisfy 


Bernoulli's Method for Complex Zeros 89 


in other words, wy and w, may be determined as solu- 


tions of the quadratic equation 
2 
x - rx+s20. (2) 


Our program generates the values mr and Sn! but 
it does not check for convergence because of lack of 
programming space. For the same reason, the starting 


x. must be supplied by the user. 


values x x 2" %3 


oO’ 


3. Flow Diagram 


Again, there is the difficulty that all available 
storage is taken up by the coefficients ay and by the 
currently needed xn° Auxiliary quantities, if they 
cannot be stored in the stack, are recalculated each 
time they are needed. The x are overwritten with Kael 


as soon as they are no longer required. 
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Display r 


Display s 
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4. Storage and Program 


00 25 = 
~ 01 RCL 3 26 - 
02 RCL 7 27 RCL 5 
03 * 28 RCL 6 
04 RCL 29 > 
05 RCL 30 RCL 6 
06 STO 7 31 RCL 7 
07 * 32 7 
08 + 33 - 
09 RCL 34 7 
10 RCL 5 35 RCL 5 
ll STO 36 RCL 6 
12 * 37 7 
13 + 38 * 
14 RCL 0 39 LAST X 
15 RCL 40 x<y 
16 STO 5 41 + 
17 * 42 PAUSE 
18 + 43 LAST X 
19 CHS 44 RCL 6 
20 STO 4 45 RCL 7 
21 RCL 5 46 > 
22 + 47 * 
23 RCL 5 48 R/S 
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5. Operating Instructions 


Load the program; switch to RUN. Select the mode of 
displaying numbers (e.g., SCI 8). Load coefficients 


ay of given polynomial into Ry (k = l, ... , 4), 
being sure to use notation (1). Load irrational star- 
ting values (such as VT, e", Y2, etc.) into Rye ae cs 
R,- Press 

PRGM 


R/S 


to start the computation. The calculator will alter- 
nate between a brief display of T and an indetermi- 
nate display of Sn: After each display of Sat 


R/S 


must be pressed to start new round of computation. 
Convergence of the sequences {ro} and {s)} must be 
tested by eye. 

If indeterminate displays of both ry and Ss, are 
desired, instruction 42 must be changed to R/S. If 
brief displays of both Th and sy are desired, instruc- 
tion 48 must be changed to PAUSE. 

An error halt may occur if an x, Or ane becomes 
accidentally zero, or if there is overflow in xy° In 
the first case, the calculation should be repeated, 
using different starting values. In case of an over- 


flow, the exponents in all x should be scaled down 
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by the same number. 


6. Examples and Timing 


Let p(x) = x4 + 3x? + 10x? = Lix? + 7. We first 


use the display mode FIX 2 (for easier reading) 
and set both instructions 42 and 48 to R/S. Using 


the starting values 


x5 ¥ R, 
x = /T > R 
Pe = e’" ” a 
a 6 
Xo := tane > Ry 
we get 
r Ss 
n n 
15.18 4.04 
- 7.05 - 35.37 
- 3.43 11.32 
- 3.85 13.11 
- 4.00 13.23 


Because convergence appears to be slow, we change 
both instructions 42 and 48 to PAUSE and let the 
computation go on. After a while, the sequences 
{x} and {s_} will have converged to within 10-2. 
We change display to FIX 5 and continue. After 


convergence to within 10> has been achieved, we 
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change display to FIX 8. After 


oe 3.97856125 , s, = 13.36969365 
has been reached, values do not change anymore. 


Dominant zeros of polynomial thus are solutions 
of 


x? + 3.97856125 x + 13.36969365 


HT] 
oO 
. 


that is, 


Wy - 1.98928063 + 3.06797266 i 
We 1.98928063 - 3.06797266 i 


el pth ig ge Oa a 
p(x) s= x + 5x + ex + Ex + 5 


After many iterations and about 18 min computing 
time, we get 
r. = 0.2756646 , s. = 0.4788911 
n n 
Thus p has the approximate quadratic factor 
2 


x - 0.2756646 x + 0.4788911 


with the zeros 


Wy 5 = 0-1378323 + 0.6781544 i 
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QUOTIENT-DIFFERENCE ALGORITHM 


1. Purpose 


To compute the zeros of a real polynomial of degree 4, 


(1) 


p(x) =a au + a,x +a 


2. Method 


The progressive version of the qd algorithm (see ENA, 
§ 8.5; ACCA I, § 7.6). For the present purpose this 
may be described as follows: One computes a sequence 


of arrays of numbers 


(n) (n) (n) (n) (n) (n) (n) 
a, 7 4, rey 1 Wo r &5 1 G3 r @y 1, 


by the initial conditions 


0 0 0 
g() 1 (0) 2 gl) 2 gl 205 (2) 


a 
pir ee ls 


o 


k= 1, 2, 3, (3) 
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and by the continuation rules 


(ntl) |_ _(n) (n) _ {(n) a 
q = qy. + ey. ee) , k=l, 2, 3, 4, (4) 
where always e(n) = e,™) := 0, and 
q intl) 
(n+l) ~_ .(n) k+l eS 
e := ey i q(nth) , k= 1, 2, 3. (5) 


k 


(Note: These formulas are most easily remembered as 


“rhombus rules": 


The system of indexing used above is not identical 
with the indexing used in the sources cited.) 


Let the zeros of p be denoted by Wy und numbered 
such that 


Then under certain conditions - for instance if the 
w, are all positive and distinct - all arrays an 
exist, and 
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lim Fy i We 7 
n--o 
k= 1, 2, 3, 4. Convergence of a to wy takes place 
: n 
with an error O(e€ ), where e« := max(|wy.i/w, |, 
Iw /wy il). By means of auxiliary computations it is 


also possible to use the arrays an to compute pairs 


of complex conjugate zeros (see ACCA I, § 7.9). 


3. Flow Diagram 


Initialize 
qy and e, by (2) and (3) 


Compute new q, by (4) 


Compute new e, by (5) 


i 


Display n :=n+1 
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Storage and Program 


Ro By By 
0 ay a5 
B 9 bart 
00 
Ol RCL 4 
02 RCL 3 
03 eu 
04 STO 6 
05 RCL 3 
06 RCL 2 
07 Z 
08 sto 4 
09 RCL 1 
10 smTo+z2 
11 RCL O 
12 CHS 
13 STO+1 
14s CL 
15 STO 0 
16 S70 3 
17. §T0 5 
18 sto 7 
+19 RCL 2 
20 STO+1 
21 RCL 4 
22 RCL 2 
23 ee 
24 sTOo+3 


RA 
w 


Quotient-Dif ference Algorithm 


5. Operating Instructions 


Load the program, and move the operating switch to 
RUN. Select the mode of displaying numbers, for in- 


stance by pressing 
FIX 8 


(adequate here because of the convergence of the q's}. 


Load coefficient ay into R. (k = 0, 1, ..- , 4). Be 


sure to define coefficients as in (1). Pressing 

R/S 
will start computation. The calculator will pause 
briefly after each array is computed, displaying its 
index n. To inspect array, press 

R/S 
during pause. The elements of the array computed last 
are then shown by pressing RCL k (k = 1, .-- , 7)- 
Pressing 


R/S 


will continue the computation. 
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p(x) = x! 


- ies" 


polynomial of order 4. The initial array 


Examples and Timing 


+ J2x" 


by the program) is 


72 
16, - 16" 0, 


The program produces the following values q 


qi 


9.39729639 
9.39507245 
9.39507092 


Time per array (including pause) 


4 
x 


p(x) 


= 10x? 


9 
7 


qi”? 


4.53453197 
4.53661877 
4.53662030 


Values produced: 


qs” 


4.0759228 
4.0040290 
4.0002261 
4.0000127 
4.0000007 


a 


2.9543419 
2.9965044 
2.9997832 
2.9999874 
2.9999993 


+ 35x 


- 96x + 24, 


(n) 
93 
1.74562397 


1.74576109 
1.74576110 


- 50x + 24 
(x - 4) (x - 3) (x - 2) (x - 1) 


ay 
.9711691 
.9994680 
.9999908 
.9999999 
.0000000 


NF PF Be 


approx. 


the Laguerre 


(set up 


(n) . 
ko? 


(n) 
q4 


0.32254767 


0.32254769 
0.32254769 


2.75 sec. 


ai”? 


0 .9985662 
0.9999986 
1.0000000 
1.0000000 
1.0000900 
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ROUTH ALGORITHM 


1. Purpose 


To determine whether a given real polynomial 


2 
p(x) := a, + a,x + ajx + ... + avx 


Wa 


of degree n 7 is stable, that is, whether ali its 


zeros have negative real parts. 


2. Method 


The Routh algorithm (ACCA I, § 6.7; ACCA II, § i-.7°. 
For the purpose of this program, this algorithm may te 


described as follows. Let 


(0) _ 
by := a 

(li 
BP Fae: 2 


k=0, 1, 2, 3. (If n< 7, a See Ba, te OT FOr 


i=l, 2, ..., ndo 
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pi) 
2 6 (eee (Gea. ie a. 
940° Ga) k dy By eR es 
0 


The polynomial p is stable if and only if the numbers 


qy" Ios ae qn all exist and are positive. 


3. Flow Diagram 


ae we pldeh) og oo fh) 
We set q := qys bh := by F by := by . To save 
storage space, quantities by are overwritten with by 


as soon as they are no longer needed. 
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4. Storage and Program 


13 - 38 
14 STO 1 39 
15 R + 40 
16 RCL 4 41 
17 i 42 
18 RCL 5 43 
19 STO 4 44 
20 : 45 
21 STO 3 46 
22 R + 47 
23 RCL 6 48 
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5. Operating Instructions 


After loading the program, switch to RUN. Load coeffi- 
cients a; of given polynomial into the storage regi- 
ster R, (i = 0, 1, ... , 7). If degree n is less than 
7, the remaining storage registers are to be filled 


with 0. Press 


PRGM 
R/S 


The computer will calculate and stop by displaying qy- 


Pressing 
R/S 


again will result in displaying G5" and so forth. The 
polynomial is stable only if the first n displayed qy 
are positive. The polynomial is unstable if some q; 

< 0 for i s n, or if display "Error" indicates divi- 
sion by zero. 


6. Examples and Timing 


p(x) = 15 + 22x + Lx? 6x? Hx, The arse 
four q,; are as follows: 
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0.79 
0.19 
0.08 


They are positive; the polynomial thus is stable. 


Indeed, its zeros are z = -l + iv2 and z = -2 +i. 


p(x) = 1l1+x + re + = + = + a + ae + x. The 


Following output is generated: 


1.00 
0.00 


Error 
Already dy = 0, indicating that p is unstable. 


SU Sub aieae eG BORN a ee eee. 


The following q,'s are obtained: 


0.17 

- 0.33 
- 24.50 
25.87 

- 0.37 
- 0.11 
0.11 
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Not all q, are positive, indicating that p is un- 
stable. Because all q exist, the fact that four 
qy are negative permits the conclusion that p has 
precisely four zeros with positive real part 
(ACCA I, § 6.7). 

Total execution time for each of these exam- 


ples is less than 20 sec. 
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SCHUR-COHN ALGORITHM I 


l. Purpose 


To determine whether or not all zeros Zs Zor Zar zy, 


of a given polynomial of degree 4, 


_ 4 3 2 
p(x) = agx + ax” + ayx” + a,x + a, (1) 
satisfy |z.] < 1. 
2. Method 


The Schur-Cohn algorithm (see ACCA I, § 6.8), in 
slightly modified form for purposes of efficient com- 


putation. For a polynomial of degree n we would form 


a triangular array of numbers a (m= 0, 1, .-. , n; 
k = 0, 1, ... , n-m) as follows: Let 
ae v= a, k = 0, l, ... , ne 


Then for m= 0, l, ... , n-l do 
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a (™) 
4 — hom 
m+1 an 
jae) ne a ~ ada as ,k = 0, 1, ... , n-m-l. 


All zeros of p satisfy |z, | < 1 if and only if the Gn 


exist and satisfy 


3. Flow Diagram 


The program implements the algorithm by means of a 
k (k = 0, l, ... , 


n) which in case of a polynomial of degree n would be 


onedimensional array of quantities b 


defined as follows: 


by am lk Sy ,™ 
v= a (™) 2 
by := ayem ¢ KE ML es yp Os 


In terms of the be the algorithm performs as follows: 


For m= 0, l, ... , n-l do’ 


eS ty by t= by - qb; (2) 
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by := by - qb. pk Sly 2) wee op Oy (3a) 

by := by qb (k-m) , k = mil, 7 n-l; (3b) 

by = bred , k=l, 2, ,n (4) 
To carry out (3), the last n-1l elements of the se- 
quence bi, bo. mais: op bon! bu Senet Je bo where bo 


occurs m-l times, are stored in the stack of the cal- 
culator. On the HP-25, this is possible for n = 4. The 
program is redundant inasmuch by is calculated m+l 


times in the m-th cycle. 
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Calculate by 


by = 


by (3), 


Pri 
at 


Yes 


Show n, Stop 


Schur-Cohn Algorithm I 


4. Storage and Program 


24 


GTO 15 


lll 
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5. Operating Instructions 


Load the program. Move the operating switch to RUN. 
Select the mode of displaying numbers. (Because it 
only matters whether lq, | <lor la, | 21, FIX 2 is 
adequate except in doubtful cases.) Load coefficient 


ay of given polynomial (1) into Ry (K-05. oD scans. ~ 5) 24). 
When 

PRGM 

R/S 


is pressed, the calculator will carry out the algo- 
rithm, briefly displaying each qy (k = 1, 2, 3, 4), 
and stop by showing n = 4. The zeros of p all lie in- 
side the unit circle if and only if all q, exist (no 
error halt due to division by 0), and if they satisfy 
la,| < 1. 

If it is desired to determine the number of zeros 
satisfying [2 | > 1, the quantities q, should be re- 
corded and used as input for the program "Schur-Cohn 
Algorithm II". Instruction 08 should be changed to R/S 
in this case, and R/S should be pressed after display 


of each q;- 


6. Examples and Timing 


p(x) = 5x4 + 4x2 + 3x? + 2x 41 
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The algorithm yields, using FIX 1 display for 


quick reading, 


Computing time 14 sec. All lay | < 1, thus all 
zeros Of p satisfy |z, | < 1. This could also have 


been inferred from the Gauss¢Lucas theorem (ACCA 
a x= 
dx x-l 
Enestrém-Kakeya theorem (ACCA I, § 6.4), because 


I, § 6.5), since p(x) = , or from the 


the coefficients a, are positive and from a mono- 


tonically decreasing sequence. 


p(x) = re + 44° + 5x" + 3x +1. 


Using FIX 1 display the program obtains 


The value 1.0 is doubtful. We repeat the calcu- 
lation with FIX 8 display, obtaining 


0.25000000 0.33333333 0.50000000 1.00000000 


Not all la, | < 1, hence not all zeros satisfy 


lz, | < 1. Indeed, p(-l) = 0. 
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SCHUR-COHN ALGORITHM II 


1. Purpose 


To determine the number of zeros Z5 of a polynomial of 


degree 4, 


4 2 
p(x) = a,x + a,x + ayX + a,x + ays (1) 


that satisfy 12, | > 1, that is, that lie outside the 


unit circle. 


2. Method 


The Schur~Cohn algorithm (see ACCA I, § 6.8). Fora 
polynomial of degree n this runs as follows: If all 
numbers q; determined by the Schur-Cohn algorithm I 
exist and satisfy Iq, | # 1, the number m of zeros out- 
side the unit circle equals me where mo is determined 


recursively as follows: 
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. I 
[nr ff dyyry! <2 
m= Bae dy Sy acah ie 2) 
Lk - Mice PE Gy aejel 2 
If some la, | = 1, the number m cannot be determined by 


this algorithm. Our program then displays m = 99. 


3. Flow Diagram 


This is a straightforward realization of the formulas 
(2). Because the qy under examination always has to be 
at the same place, the remaining q, are shifted up one 
place at the end of each cycle. The program implements 


the case n = 4, 
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Display "99" 
Stop 


Display m 
Stop 
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4. Storage and Program 


Ro Ry R, R, R, Re Re Ro 
qy qo q3 q, k m 
00 25 STO 2 
ol CLX 26 GTO 04 
02 STO 6 > 27 9 
03 STO 7 28 9 

> 04 1 29 GTO 00 
05 STO+6 + 30 RCL 7 
06 RCL 4 31 GTO 00 
07 ABS 32 
08 x <y 33 
09 GTO 16 34 
10 x = y 35 
ll GTO 27 36 
12 RCL 6 37 
13 RCL 7 38 
14 = 39 
15 STO 40 

+ 16 RCL 6 41 
17 4 42 
18 x = y 43 
19 GTO 30 44 
20 RCL 3 45 
21 STO 4 46 
22 RCL 2 47 
23 STO 3 48 
24 RCL 1 49 
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5. Operating Instructions 


Load the program and move operating switch to RUN. 


Press 
FIX 0 


for integer display of number of zeros. Load qs into 
R., i=1, 2, 3, 4. When 


PRGM 
R/S 


is pressed, the calculator after a short time will 
display the number m of zeros outside the unit circle. 
If some la, | = 1 and hence the algorithm fails, the 
integer 99 will be displayed. 


6. Examples and Timing 


With the data of example of the preceding 


algorithm we get 


Indeed, we already know that this polynomial has 
no zeros outside the unit circle. Computing time 


about 4 sec. 
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p(x) = 27x) - 18%? + 52x? + 3x - 15. 
The program “Schur-Cohn Algorithm I" yields the 


values 
qa = -0.6, -0.4, 4.7, 0.2 
Using these data, the present algorithm yields 


m = 2, indicating that precisely two zeros of p 


satisfy lz, | > 1. 


Part 4 
POWER SERIES 
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RECIPROCAL POWER SERIES 


1. Purpose 


Given a formal power series with leading coefficient 1l, 


P=1+ x +a ae + 
= ay 2 genes 


to compute the coefficients bo of the reciprocal power 


series 


-1l 2 
P = 1+ b)x + bx 5 a 


satisfying p+ P=1. 


2. Method 


This is one instance where the very limitations of the 


pocket calculator make it necessary to invent a new 
To ee 


ee 


and unorthodox algorithm to solve the problem. Ordina- 


rily one would use the fact that 


2 2 = 
(l + a,x + a,x” + wee) (1 + byx + box +...) = 1, 
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n lvn-1 2°>n-2 n ‘ 
to obtain the recurrence relation by = l, 
Dig eo Ae acy agPa VS anPo , 
n=1, 2, ... . However, to compute bo by this method 


requires computing and storing all be where k <n, 
which exceeds the storage capacity of the pocket cal- 
culator except for trivial values of n. The same is 
true for "fast" algorithms for computing the recipro- 
cal series using Newton's method and the Fast Fourier 
Transform (see Aho, Hopcroft, and Ullmann, The design 


and analysis of computer algorithms, Addison-Wesley, 
1974). 


Instead, we use the following approach which, how- 


ever slow the resulting algorithm, works on the HP-25. 
Writing 


so that P = 1 + Q, we obviously have 


Pp =1-Q+4 9? - Q? + ... 


For any n 1, the coefficient bo, thus is the sum of 
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all coefficients of the power x” in all powers o* 
where 1 = k n. In a fixed power o*, the coefficient 


n ‘ 
of x arises as the sum of all products 


h 
(-l) a a bag a ° 
aie Ky 
where k, =1l, i=l, 2, . , nh, and 
ky + ky t+ ..e F ky =n 


Y(-ay dees (rap ds (1) 


where the sum comprises all systems of indices ky, koe 


.++ , kp such that h 24, 


liv 


(2) 


(No upper limit on h needs to be imposed.) 

At first sight, it seems difficult to programa 
complicated formula such as (1) for a pocket calcula- 
tor, due to the many systems of indices (Ky rkoreeerky) 
that have to be kept track of. However, the following 


device, suggested to the author by E. Specker, makes 
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the programming possible. The systems of indices k = 
(kK) okos--- ky) satisfying (2) are put into a one-to- 


one correspondence with the integers m satisfying 


nod s Im < a in the following way: Let 
= 2. n-l 
m= £5 + 2c) + 2 Ee» + oes F 2 fo. (3) 
be the binary representation of m (€. = 0 or 1, j = 0, 
1, ..- , n-2; en-1 = 1). Por each j such that oe = 1, 


we put a mark at the point x = j of the real line. In 
addition, we put a mark at x = 0. These marks divide 
the interval [0, n-1] into subintervals. The indices 


k, of the system (ki +k sky) corresponding to m are 


gree 
the lengths of these subintervals. 
Accordingly, the coefficients bo May be calculated 


by the formula 


where the products qi are defined by means of the bi- 


nary expansion of m, 


here k k, are the gaps between the non- 


1’ Kae oe 
zero ey in (3). 
The binary expansion of m is calculated by the 


following algorithm. Let Pp t= 


Reciprocal Power Series 127 


Ps. 
Be BS CIN 4S i De eet te 
3 2 
Then 
Pp. 
1, if FRAC 5240, 


0, otherwise. 


(m) 
(0) 


thus is calculated re- 


:= 1, 


Given m, the product I := IN 


cursively as follows: Let I 


163) e.=0, 
(341) 8 


I 
(3) ; 
eeuct er 


here k := j - j' where j' := max { i [ i<j, e, #0 }. 
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3. Flow Diagram 
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4. Storage and Program 
R R R R R R R R 


0 1 2 3 4 5 6 7 
or ree p n [0 
m 2 x 
00 25 
01 Ll 26 
02 STO 6 27 
03 STO-1 28 
04 STO-2 29 
05 RCL 1 30 
06 STO 5 31 

> O07 CLX 32 
08 STO 3 33 

» 09 l 34 
10 STO+3 35 
ll RCL 5 36 
12 2 37 
13 a 38 
14 INT 39 CHS 
15 STO 5 40 STO*6 
16 LAST x 41 RCL 5 
17 FRAC 42 x #0 
18 x =0 43 GTO 07 
19 GTO 09 44 RCL 6 
20 45 STO+7 
21 46 RCL 2 
22 47 x #0 
23 48 GTO 01 


24 49 RCL 7 
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Locations 20 + 38 are reserved for the program for the 
computation of ay. This program should assume k in R3i 
at the end of program execution, ay must be in X re- 
gister. The registers Ry and R, are available for 
auxiliary storage. 

5. Operating Instructions 

Load the program, including the program to compute aye 


Unused locations must be filled up by NOP instructions. 
Move the operating switch to RUN. Select the mode of 
displaying numbers. If coefficient bo is desired, load 


data as follows: 


2 into 
2°71 into 
0 into R 


Here the numbers 2n and 2n-l must be genuine integers; 
the calculator will not deal correctly with a 2” that 


is computed by the y* instruction. Press 


PRGM 
R/S 


to start computation. After some time (see examples) 
the calculator will display bie Caution: bo is calcu- 


lated as a sum of 2” terms which may differ in sign 
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and size. Thus already for moderate n, bo may be af- 


fected by rounding error. 


6. Examples 


To calculate the coefficients bo in the expansion 


= 1+ b,x + bax’ te casera’ 3 


a 
Log(l - x) 1 


(The bo are required for certain methods of nume- 
rical integration; see ENA, § 13.4.) The fore- 


going series is the reciprocal of 


Po:= hog ()_~_x) ee er ee ae 


are generated by means of the program 


20 RCL 3 24 

21 1 : 

22 ‘ : NOP 
23 1/x 38 


Results: 
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Bb (calculated) bo (exact value) 
-0.500000000 = 7 = -0.500000000 
-0.083333333 a GS = -0.833333333 
-0.041666667 - ai = -0.041666667 
-0.026388889 - 2 = -0.026388889 
720 : 
3 L 
-0.018750000 - Teo = 0.018750000 
863 
-0.014269180 - 60480 ~ -0.014269180 


Computing 


time 


109 


252 


sec 


sec 


sec 


sec 


sec 


sec 


All calculated values are accurate in all digits 


given. 


To compute the Bernoulli numbers By 
x = 5 By k 
ar Pie — x 
e -l1 k=0 * 
(see ACCA I, p. 13). - Setting 
x 
e -l_ 2 
Sager ee 1+ a,x + a5x Be yea 


the program will calculate the numbers 


defined by 


Reciprocal Power Series 


The coefficients 


La 2 
k (k+l)! 


may be calculated by the following program: 


20 1 28 GTO 23 

21 STO+3 29 RCL 3 

22 RCL 3 30 1/x 

23 1 31 

24 - 

25 x = 0 : NOP 

26 GTO 29 ; 

27 STO*3 38 
Results: 


bo (calculated) bo (exact values) Computing 


time 
-0.500000000 - = = -0.500000000 3 
0.083333333 T= 0-083333333 8 
ee ee 0 20 
~0.001388889 - 3 = -0.001388889 54 
goa too 0 129 

? 


0.000033069 303 


0.000033069 


sec 


sec 


sec 


sec 


sec 


sec 
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POWER OF POWER SERIES 


1. Purpose 
Given a formal power series with leading coefficient l, 


2 
P=1+ a,x + a,x + asx + eee, (1) 


and given a real number a, to compute the coefficients 


ate of the a-th power of P, 


a 2 oy 
(1 + ayx + ax + ...) 


ae] 
i] 


= (a) (a) 2 
= 1+ ay x + as Mo FS oes 


2. Method 


By conventional wisdom, the purpose would be served by 
means of a recurrence relation derived from the fact 
that the series R := po satisfies the formal differen- 


tial equation 
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(J. C. P. Miller formula; see ACCA I, Theorem 1.6c). 
This approach is not feasible for pocket calculators 
due to storage limitations. Instead, we compute the 
series Pp“ directly from its definition as composition 


of the binomial series 


with the series 


This yields 


hence 


)(-ay ) + (ray ), (2) 


where the sum comprises all systems of indices (kik, 

++ 7K) satisfying the conditions (2) of the preceding 
algorithm. These systems are coded in terms of the bi- 
nary representation of the integers m satisfying gn-1 


< n 
=m < 2°, and we obtain 
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li, (3) 


where qa has the same meaning as in the preceding al- 
gorithm, and where h is the number of the non-zero di- 


gits in the binary representation of m. 
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3. Flow Diagram 


Power Series 


- Storage and Program 


Ro a Re 

9m gn-l 
m g 
00 

> 01 l 
02 STO 6 
03 sTo-1 
04 STO-2 
05 RCL l 
06 STO 5 
07 CLX 
08 STO 4 

+ 09 CLX 
10 STO 3 

> ll 1 
12 STO+3 
13 RCL 5 
14 2 
15 + 
16 INT 
17 STO 5 
18 LAST x 
19 FRAC 
20 x = 0 
21 GTO 11 
22 
23 
24 
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Locations 22 : 31 are available for the program to 


compute a This program should assume k in Ry at the 


k* 


end of program execution, a, must be in X register. 


k 
(If it is more convenient to compute l/a,, this may be 
done; instruction 32 must then be changed to STO:6.) 

No registers are available for auxiliary storage, but 


stack may be used freely. 


5. Operating Instructions 


Load the program, including the program to compute ay 
Unused program locations are to be filled with NOP in- 
structions. Move the operating switch to RUN. Select 


the mode of displaying numbers. Load exponent 


a into Ro 


and, if coefficient af) is desired, 


2 into R 
n-1 : 

2 into R 

0 into R 


(The powers of 2 should be loaded as integers, without 


using the y* instruction.) Press 


PRGM 
R/S 
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to start computation. The calculator will stop by dis- 


. a “ F : : 
playing a! ns Caution: Computing time increases expo- 


nentially with n. For large n, a0 may be contamina- 


ted by rounding error. 


6. Examples and Timing 


Power series reversion. Let y be given by a 


power series in x of the form 


Then the expansion of x in terms of y is 
2 3 
x= y+ Coy + Cay + eee 


where, by the Lagrange-Biirmann Theorem (see ACCA 
I, § 1.9) 


-n = 
Here res R denotes the coefficient of x : in 

‘ -n 
the expansion of R.. Because R = xP, where P has 


the form (1), this coefficient in the foregoing 
(-n) 


i e 1 
notation equals an-l 


There follows 


a7?) 


=l ‘“ 
Co =n And 7 2p Sy ees 


n 
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As an example, consider 


Ree -l, 
where 
Ae ee ott 
k (k+1)! 
Because 
o n-l 
a Ul exteg aia ae i ge 
n 
n=1 
we have in this case 
n-1 
_ (-1) (-n) _ ,_,,n-1 
o, ~ n e-Bay. 


The following program is used to compute the aye 


22 1 28 GTO 31 
23 STO+3 29 STO*3 

24 RCL 3 30 GTO 25 
25 1 31 RCL 3 

26 7 32 STO=6 

27 x = 0 


Being ignorant of the exact answer, the main pro- 


gram then computes the following values: 
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n ace Computing time 
2 -1.000000000 4 sec 
3 1.000000000 10 sec 
4 -1.000000000 23 sec 
5 0.999999999 58 sec 
6 -1.000000000 137 sec 
7 0.999999998 5.5 min 
8 -1.000000006 12 min 
9 1.000000008 

10 -0.999999858 

ll 0.999999631 


Stirling's formula. The coefficients c, in the 
asymptotic expansion of the [ function, 


x-y2 c “3 


T(x) % /2n x e*~ (1+ + fee a ladip OE Py 


can be defined as 


(see ACCA II, § 11.6). The program calculates the 
following values: 


Power 


(-4-3) 
Pog q 


-083333334 0 .083333334 


-001157408 0.003472223 


-000178759 -0.002681385 


-000002234 -0.000234570 - 


of Power Series 143 


exact value 


se 
D- 0.083333333 
—— 0.003472222 
288 020 * 
139° _ | 
- Siea0 = 0.002681327 
571 


3488320 -0.000229472 
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EXPONENTIATION OF POWER SERIES 


1. Purpose 


Given a formal power series with constant coefficient 


zero, 


to calculate the coefficients Ce of the series 


Q_ 2 
e~ =l + cx + Cox Bee we 


2. Method 


To use the differential equation satisfied by the se- 
ries S := e® is not feasible due to the lack of stora-~- 


ge facilities. However, from the fact that 


we have, using the notation of the two preceding pro- 


grams, 
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where the sum comprises all systems of indices (ky ky, 
veerky) Satisfying the conditions (2) of the program 
"Reciprocal Power Series". As before, these systems 
are coded in terms of the binary representation of the 
integers m satisfying ane sm< rae and we thus ob- 


tain 


where Te has the same meaning as in the two preceding 


programs. 
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3. Flow Diagram 
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4. Storage and Program 


Ry Ry R, R, R, Re Re Ra 

on gt! og h P R [0] 
m 2 z 
00 25 

> Ol 1 26 
02 STO 6 27 
03 STO 4 28 
04  sTo-1 29 
05 sTo-2 30 
06 RCL 1 31 
07 sTO 5 32 

> 08  CLX 33 
09 STO 3 34 

> 10 1 35 
11 = STO+3 36 STO*6 
12 RCL 5 37 RCL 4 
13 2 38 STO+6 
14 : 39 1 
15 INT 40 STO+4 
16 sTO 5 41 RCL 5 
17. “LAST x 42 x #0 
18 ~~ FRAC 43 GTO 08 
19 x = 0 44 RCL 6 
20 GTO 10 45 STO+7 
21 46 RCL 2 
22 47 x #0 
23 48 GTO Ol 
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Locations 21 + 35 are to be used for the program to 
compute ay. This program should assume k in R53 at the 
end of program execution, ay must be in X register. 
(If it is more convenient to compute l/ay, this may be 
done; instruction 36 must then be changed to STO+6.) 
Register Ro is available for auxiliary storage. Unused 


locations are to be filled with NOP instructions. 


5. Operating Instructions 


Load the program, including the program to compute aye 
Move the operating switch to RUN. Select the mode of 


displaying numbers. If coefficient ch is desired, load 


2 into R 
n-l 1 

2 into R 

0 into RL 


(The powers of 2 should be loaded as integers, without 


using the y* instruction.) Press 


PRGM 
R/S 


to start computation. The calculator will stop by dis- 
playing oe (Computing time increases exponentially 
with n.) Caution: For large n, if some ay are negative, 
cancellation may cause c, to be contaminated by roun- 


ding error. 
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6. Examples and Timing 


Q) 


Q = Log(l - x) 


We should obtain 


that is, cy 


to compute a 


Results: 


wm f& WwW Ne 


k 


« lk 
=- J =x 
uo 
2 =l-x, 
-l, e..= 0 for k > l. 


is as follows: 


21 RCL 3 
22 CHS 
23 GTO 36 
24 
. NOP 
35 
36 STO+6 

c 

n 


-1.000000000 
0.000000000 
0.000000000 

1 * 19720 
0 .000000000 


Computing time 


3 
6 
15 
38 
90 


sec 


sec 


sec 


sec 


sec 
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The program 


150 Power Series 


0 .000000000 
1 * 19°19 


3.0 min 


8.1 min 


The Bell numbers b_ are defined by the expansion 


For 


we may calculate a 


21 RCL 3 28 
22 1 29 
23 - 30 
24 x = 0 31 
25 GTO 28 
26 STO*3 
27 GTO 22 36 
Results: 
n c b =nic 
n n n 
1 1.000000000 1 
1.000000000 2 
3 0.833333333 5 


i x 
k: 


k 


k by the program 


RCL 3 
STO+6 
GTO 37 


NOP 


Computing time 


3 sec 
8 sec 


18 sec 


wow on Dn Ww & 


Exponentiation of Power Series 


0.625000000 
0.433333333 
0.281944445 
0.174007996 
0.102678571 
0.058275463 


15 

52 
203 
877 
4140 
21147 


46 
110 
260 

10 

23 

51 


sec 
sec 
sec 
min 
min 


min 
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Part 5 
INTEGRATION 
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NUMERICAL INTEGRATION WITH STEP REFINEMENT 


1. Purpose 


To evaluate the integral 


numerically for an arbitrary function f that can be 


evaluated by a program requiring no more than 16 in~ 
b 

structions. The general integral f is reduced to the 
a 


above by the substitution x' = x a. 


2. Method 


We approximate I by a sequence of approximate values 
Ihe k = 0, 1, 2, ... . The approximation I, is ob- 
tained by dividing [0, b] into 3% congruent subinter- 
vals and evaluating the integral on each subinterval 
by the midpoint formula. (The trapezoidal rule is 
avoided because integrands frequently have singulari- 


ties at the endpoints of the interval. Although harm- 
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less for the existence of the integral, these may 


cause the program to fail; see example .) Writing 
_ 47k 
hy 7= 2 “b, 
we have 
274 : 
I. = hy z £((m + 5) hy) 5 
m=0 


For any continuous f, 


lim I, =I 


kro 


fff is sufficiently smooth, the convergence of the 
sequence {I,} May be sped up by the Romberg algorithm. 
To this end the program saves the five currently most 
recent values of qT, in locations in which they can be 
used as input for the Romberg acceleration program 


(see the following program) without external storage. 
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3. Flow Diagram 


compute f(x) 
s + £(x) 
x +h 


Display 


push down I 


k 


zh 


158 


4. 


Integration 


Storage and Program 


R 


0 


> 


STO 4 
PAUSE 
2 
STO+5 
RCL 1 
STO 0 


(R/S) 


oO MW BN WwW Ft WN 
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The program to compute f(x) should be in locations 33 


through 49, assuming x in R Only the stack may be 


6° 
used for temporary storage. The last instruction must 
be GTO 09. At this point, f(x) must be in the X regi- 


ster. 


5. Operating Instructions 


Load the main program. Load the program for computing 
f into locations 33 through 49; the last instruction 
must be GTO 09. Switch to RUN. Select a mode of dis- 


playing numbers, for instance 
FIX 8 


If f involves trigonometric functions with argument i: 


radians, press 


To start computation, press 


PRGM 
R/S 


The calculator will pause briefly while displaying 


each I,. If k 2 0, and if 


R/S 
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is pressed during pause, the values I will be found 


k-m 
in Ry _ in! m= 0, l, 2, 3, 4. 


6. Examples and Timing 


b=, £(x) := _ . Program to compute f: 
l+x 

33 4 

34 RCL 6 

35 x? 

36 1 

37 + 

38 + 

39 GTO 09 


The following values of I, are obtained: 


k 
(k) 

k 1 T, 

0 3.20000000 } 3.14191817 
1 3.16235294 3.14159265 
2 3.14680052 + Romberg 3.14159264 
3 3.14289473 3.14159266 
4 3.14191817 { 3.14159264 


Time required 35 sec. Subjecting the values I, 


to the Romberg algorithm produces the exact value 


nm = 3.14159265 with an error of one unit in the 
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last digit due to rounding effects. 


Values obtained: 


k 

0 0 
1 0 
2 0 
3 0 
4 0 


Ty 


-81093022 
-81936429 
-82167416 
-82226766 
-82241711 


N 


aL 


} 


b= 1, £(x) := wept tt kt . Program to compute f: 
33 1 
34 RCL 6 
35 + 
36 ln 
37 RCL 6 
38 
39 GTO 09 


-82241711 
-82246693 
-82246702 
-82246702 
-82246702 


Romberg. 


oo 00 8 


Exact value: l2 = 0.82246703. Time required to 


generate [I 


orcte et 


I 


4 


about 40 sec. 
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ROMBERG ALGORITHM 


1. Purpose 


To accelerate the convergence of a sequence {Tv} to 


its limit age where 


ae := T(2 ™h) , m=0, 1, ..-, 4, 


nd the function T(h) has an asymptotic expansion of 


the form 


2 4 6 
T(h) % ay t+ ayh” + ajho + ash +... , ho» OF. 


This problem poses itself in connection with numerical 
integration [?(h) = approximate value of integral cal- 
culated with step h]; in addition, it occurs in nume- 


vical differentiation, and in the computation of cer- 


tain functions and constants. See examples below. 


2. Method 


The Romberg algorithm (see ENA, § 12.4; ACCA II, 
§ 11.12). With the data 
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OC ea a ee ee 


we generate a triangular array of numbers ao by the 


formulas 


If m were allowed to tend to infinity, one would have 


(loc. cit.) for each fixed n 


mn 


+ 0(4.) 


Thus Ty, may be regarded a much more accurate approxi- 


mation to ay than Tyo (mathematical error bounds are 


available). To check convergence, the program is ar- 


ranged so that all values T (n = 0, l, ... , 4) are 


4n 
available at end of computation. 
For reasons of numerical stability as well as eco- 


nomy, formula (1) is evaluated as 
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3. Flow Diagram 


Because an address modification is not available on 
the HP-25, the operation (2) is always carried out at 
the same place. This is made possible by a cyclic per- 


mutation (here called "rotation") of the Th 


Rotate 


:= mt+l 


Display Tag 


Stop 
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4. Storage and Program 


Ro Ry R, R3 Ry Re Re Ro 
a, Tr T Ts rs temp n m 
00 25 RCL 1 
4 ol 1 26 sTo 0 
02  sTO 6 27.» -RCL 2 
ae ee 28 sto l 
_04 sto 7 29° RCL 3 
> 05 RCLG6 30 sTo, 2 
06 RCL 7 31 RCL 4 
07 x<y 32. «STO 3 
08 GTO 23 33. RCL 5 
4 09 4 34 STO 4 
10 RCL 6 6° 35 1 
ll y™ 36 STO+7 
12 1 a 37. RCL 7 
13 - 38 4 
14 1/x 39 x =y 
15 RCL 1 _ 40 GTO 05 
16 RCL O 8 41 1 
17 - _ 42 _ STO+6 
18 * J 43 RCL 6 
19 RCL 1 44 5 
20 + 45 xy 
21 sto 0 46 GTO 03 
__22 PAUSE (R/S) AO 47 RCL O 
§ + 23 RCL O 48 GTO 00 


24 STO 5 49 
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5. Operating Instructions 


Load the program, and switch to RUN. Choose the de- 


sired mode of displaying numbers, for instance 
SCI 8 


to get floating eight-digit display. Load data Th = 
T(2°™n) into Rae m= 0, l, ... , 4. (If the data are 
produced by the program "Numerical Integration with 


Step Refinement", they are already there.) When 


PRGM 
R/S 


is pressed, the calculator will generate the array Ta 


column by column, pausing briefly after computing each 


entry. (For an indeterminate display of each Te ’ 


n 
change instruction 22 to R/S and press R/S after each 


display.) The calculator will stop by displaying Tygr 


Ty, is also in Roi generally, the elements in the last 
row of the array, T, n! are found in the registers 

, 
Rain! n= 0, l, ... , 4. 


6. Examples and Timing 


Numerical evaluation of 


Romberg Algorithm 
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by the preceding program produces the approximate 


values 


0.001953125 
0.037544250 
0.078839094 
0.094288322 
0.098544474 


The present program accelerates the convergence 


of these values to produce the following values 


T 


4,n 


in the last row of the Romberg array (listed 


vertically here for reasons of space): 


0.098544474 
0.099963191 
0.099998199 
0.099999858 
0.100000000 


It is seen that the last value is accurate to all 


places given. Computing time about 44 sec. 


Here we compute Euler's constant, 


Y 


:= lim (1 + 


n> 


Nie 


BB aera 


+ BS - Log n) 


’ 


Integration 


directly from its definition, using Romberg acce- 


leration. If for n= l, 2, ... we let 
= i ,i ee ee 

s(n) := 1 + 5 + 3 t+... + a. + 55 Log n, 
then we evidently also have y = lim s(n). The 
factor 5 in front of is motivated by the fact 
that for suitable ay 

s(n) ®y + £ ayn >* j n+ 

k=1 


(see ACCA II, § 11.11). Thus the numbers 
k 
% i= s(2°) , k = 0, l, 2, ... , 
satisfy the conditions for applying the Romberg 


algorithm (see ACCA II, § 11.12). They are con- 
veniently generated by the recurrence relation 


Sy41 = I F 


k = 0, 1, 2, ... . It iS an easy exercise to 
write a program that computes 9; for i= 0, l, 


--. , 4 and stores it in Rj: We obtain 
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0.500000000 
0 .556852820 
0.572038973 
0.575915602 
0.576890272 


The program "Romberg Algorithm" accelerates these 


values as follows: 


0.576890272 
0.577215162 
0.577215652 
0.577215663 
0.577215665 


Again the last value is accurate to all places. 


Computing time for the acceleration about 44 sec. 


For further examples on Romberg's algorithm see 
the programs "Numerical Integration with Step 
Refinement", "Plana Summation Formula", “Log- 


Arcsine Algorithm". 
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PLANA SUMMATION FORMULA 


1. Purpose 


To evaluate the integral 


by means of numerical integration. This integral oc- 
curs, among other places, in the Plana summation for- 
mula (ACCA I, § 4.9), 


wo 


fy ae Hie ax ey SOO ae 
0 0 


Xo £(n) = an 
= e Yy _ 1 


fe) 


Nie 


n 
where 
gfy) := - 2 Im f(iy) . 


This formula holds, roughly, for functions f that are 


2uIlv1 gor 


analytic for Re z 2 0 and grow less than e 
y := Im z-+ to . It may often be used to sum slowly 
converging series. The program works for any function 
g whose evaluation requires no more than 18 instruc- 


tions. 
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2. Method 


The integral I is approximated by the sums 


(@) g(yy) 
In = ¢£ an , m= 0, l, 2, ... , 
k=0 Y 
e - il 
-m h : , 
where h = ho - 2 *: Yx = 5 + kh. The summation is 
any 
stopped when e Ride 1 > l/e. The constants hy and l/e 


are input. The program displays each Th after it has 
been computed. The convergence of the sequence {I} 
may be accelerated by the Romberg algorithm (see the 


preceding algorithm). 
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3. Flow Diagram 
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4. Storage and Program 


Ro Bal BD Ry a Re Re By 
h l/e g z y 
00 25 1/x 
> Ol RCL O 26 RCL 6 
02 2 27 2 
03 J 28 T 
04 sto 2 29 * 
05 STO 6 30 * 
06 0 31 e* 
07 $T0 5 32 1 
> 08 RCL 6 33 - 
09 3 34. RCL 1 
10 » ?P 35 x < y 
11 In 36 GTO 43 
12 > P 97 R + 
13 x? 38 + 
14 xy 39 sTo+s 
15 2 40 RCL O 
16 * 41 STO+6 
17. —- RCL 6 42 GTO 08 
18 3 43. RCL O 
19 +P 44 RCL 5 
20 R + 45 * 
21 + 46 R/S 
22 sin 47 RCL 2 
23 + 48 STO 0 
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Note: The function g is evaluated by instructions 08 
through 25, presuming its argument y in Re: gly) 

should be in X register before executing instruction 
26. In above program, g is taken from example be- 


low. 


5. Operating Instructions 


After loading the program, switch to RUN. Indicate the 


mode of displaying numbers desired, for instance, 
SCI 8 
If function g involves trigonometric functions with 


arguments in radians (such as in examples [2] and 
below), press 


Load hy (such as 0.5) and press 
STO 0 
10 
Load 1/e (such as 10 = 10) and press 


STO 1 


When 
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PGRM 
R/S 


is pressed, the calculator will compute and display 


I,. Pressing R/S again yields I I 


0 1 ae ene 


The numbers I are not stored. To subject 


I 
Oe A 
them to the Romberg algorithm, they should be copied 


and used as input for the program Romberg Algorithm. 


6. Examples and Timing 


Euler's constant y (see preceding algorithm) may 


also be represented as 


ao 
veges 4, = ey 
0 l+tyie - 1 


(ACCA I, p. 275). The function 


2 
g(y) := —, 
l+t+y 
is generated by 
08 RCL 6 
09 2 
10 * 
ll RCL 6 


12 x 
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13 

14 + 

15 = 

16 GTO 26 


The following values Ih result by choosing ho := 
0.5, lfe := Lore: 


6.6297262 * 10 
7.4582004 * 10. 
7.6562828 * 10 
7.7052792 * 10° 
7.7174967 * 10 


Total computing time: approx. 10 min. Subjecting 
the foregoing values to the Romberg algorithm 
yields the following last line of the Romberg 
scheme: 


7.7174967 7.7215692 7.7215663 7.7215664 7.7215664 


(all * 1077). We conclude that I = 7.7215664 * 


10°) and thus 
y = 0.577215664 , 


which is correct to within a rounding error of 


1* 10°. 
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Application of the Plana formula to f(z) := 
(l + z) * (principal value) shows that the Rie- 


mann zeta function (ACCA II, § 10.8), 


G(s) := lS 
n=l on 


may be represented as follows: 


where 


.. 2 sin(s Arctg y). 
(1 + gore? 


This representation holds for all s # 1. We load 


s into R, and use the following program to com- 


3 
pute g: 
08 RCL 6 
09 tan? 
10 RCL 3 
1l * 
12 sin 
13 2 
14 i 
15 RCL 6 


16 x 
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17 1 
18 + 
19s RCL 3 
20 y* 
21 vx 
22 + 


For s := 1.1, hy 2= 0.5, lfe := 107° we obtain 
7.2387260 * 1077 } f g.4403697 * 10° 
8.1548623 * 10° 8.4448497 * 10-7 
8.3730175 * 107° Romberg, < 8.4448464 + iors 
8.4269295 * 107 8.4448464 * 102 
8.4403697 * 1077 8.4448464 * 107° 


There follows ¢(1.1) = 10.584448464. The program 
may be tested by evaluating ¢(s) at the "trivial" 


zeros s = -2, -4, 


To evaluate the very slowly converging series 


> 1 
Si= f ————,; 
n=3 n{Log n} 


, 


we use the Plana formula where 


L 


f(z) <= Paes eG Oe 
(3 + z){Log(3 + 2)} 
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This requires the evaluation of 
g(y) := - 2 Im f(iy) , 


which may be written as 


_ 2 sin(a + 28) 
g(y) =a 
rR 


where 


3 + iy=r e , Log r + ia a8 


" 
a 
o 


The program is given in the foregoing main pro- 
gram; it offers interesting applications of the 
» P operation. The final multiplication by 2 is 


omitted for lack of space. The following values 


result with ho := 0.5, lfe := 101°, 

0.8983867 * 1077 f 1.0359638 * 10° 
1.0025014 * 10° 1.0364923 * 1072 
1.0280286 * 10-2 Romberg, { 1.0364921 * ied 
1.0343783 * 1077 1.0364921 * 1077 
1.0359638 * 10 2 1.0364921 * 1077 


We conclude 


l 1 
3 (Log a) 


+ 2 * 1.0364921 * 1072 


Ni 


Log 3 


1.0690583 
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The total time required for this computation is 


about 20 min. 
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DIFFERENTIAL EQUATION OF FIRST ORDER: 
TRAPEZOIDAL METHOD 


1. Purpose 


To integrate numerically a general first order equa- 


tion 
y' = £(x,y) (1) 
subject to the initial condition 


y (Xo) = Yo - 


2. Method 


The trapezoidal rule. Choosing an integration step h, 


we determine approximate values Yn of the values y(x,) 


of the exact solution y(x) at the points X, 1 Xp + nh 
(n = 1, 2, ...) by the formula 
-y = Bletx yy.) + £x ey |. (2) 
Yn+1 n 2 n/7n n+1/7n+1 


The implicit equation for Yne1 that arises at each 
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step is solved by iteration, using Euler's method as a 
predictor (see ENA, § 14.7). Iteration is stopped when 
two successive iterates differ by less than c€. Eleven 
storage locations are available to evaluate the ex- 
pression DE (x,y) - 

The order of the method defined by (2) is 2, inde- 
pendently of the number of iterations performed. That 
is, as h + 0, the approximate values Yo at a fixed va- 
lue of x tend to the exact value y(x,) with an error 
of Oth?) . If « is small (say, e e0o"*) and the equa- 
tion (2) thus is satisfied to full machine accuracy, 
the convergence of Yn to y(x)) may be sped up by the 
Romberg algorithm (see the example). However, a large 
number of iterations may then be necessary at each 
step. The number of iterations may be reduced by gi- 
ving € larger values. For very large values of ¢« (Say, 
€ = 10?°), the corrector will be applied only once, 
and the method becomes identical with the simplified 
Runge-Kutta method, 


_h 
Yona YR = BEG ry,) + OK gyrate Oy r¥q))] (3) 
The order of the method is still 2, but Romberg's noe 


extrapolation is no longer feasible. 


3. Flow Diagram 


We write x := Xa y:= Yni y* denotes the current ar- 
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gument y in f(x,y), y* is the corrected value. 


Display x 
Display y 


A “flag" k is used to distinguish between the predic— 


tor (k = 0) and the corrector cycles (k = 1). 
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4. Storage and Program 


» 05 


Program to 


compute 


£(x,y*) 


RCL 0 


> 
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>. Operating Instructions 


Load the main program. Load the program for computing 

f£(x,y*) into program locations 05 : 15, assuming x in 
Pa 

R, and y* in RL. 

X. (As an alternate possibility, one may use locations 


This program must place f£(x,y*) into 


05 + 20 for a program that places A := BE (x,y) into 


R,.) R. is available to store a constant required for 


6° 5 
the evaluation of f. 

Move the operating switch to RUN. Select the de- 
sired mode of presenting results, for instance by 


pressing SCI 6 to get floating six-digit display. Load 


integration step h into Ro 
starting value Xo into Ry 
starting value Yo into Ry 
tolerance e« for stopping 
inner iteration aRES Ry 
When 
PRGM 
R/S 


is pressed, the calculator will start computing. After 
completing inner iteration, it will briefly display x 


and stop by displaying Yy° Pressing 


1 


R/S 
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will cause the display of Xr Yo: etc. The process may 
be repeated until the desired range of Xn has been co- 


vered. 


To check the accuracy of the results, the computa- 


h 


tion should be repeated with the steps t, qrocee un- 


til the results have converged to the desired accura- 
cy. If e€ is small, the convergence may be sped up by 
the Romberg algorithm. 


6. Example and Timing 


The exact equation of the mathematical pendulum of 
length 2, 


QR 
(g := gravitational constant = 9.81 ms?) may under 
the initial condition y'(0) = 0 be integrated to give 
y'=- jeg Ycos y - cos y 
2 0’ 


where Yo = y(0). We solve this equation under the ini- 


tial condition 
T 
Yo = y(0) = 37 1.047198 . 


Here 
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=f = ivy f/f LD Seeety = 0L5 
A = 5 f(x,y) = h 7h cos y 0.5. 


The following program is used to compute A: 


eke) RCL 7 
06 cos 
07 : 
08 5 
09 = 
10 vx 
1l RCL 5 
12 * 
13 RCL 0 
14 > 
15 CHS 
16 STO 6 
17 NOP 
18 NOP 
19 NOP 
20 NOP 


The constant / 54 is assumed in Re: Since the argument 
of cos is assumed in radians, RAD must be pressed be- 
fore starting computation. 

The following values result if £ = 1, beginning 


with h = 0.1 and successively halving the step: 
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Y 
x h=0.1 h=0.05 h=0.025 h=0.0125 
0.0 1.047198 1.047198 1.047198 1.047198 
0.1 1.005241 1.004979 1.004913 1.004896 
0.2 0.881649 0.880558 0.880284 0.880215 
0.3 0.683836 0.681290 0.680648 0.680488 
0.4 0.425723 0.421160 0.420008 0.419719 
0.5 0.128475 0.121689 0.119976 0.119546 
0.6 ~0.180863 -0.189433 -0.191595 -0.192137 


Applying Romberg to the values at x = 0.5 yields 


-128475 

-121689 0.119427 

-119976 0.119405 0.119404 

-119546 0.119402 0.119402 0.119402 


oo oOo 8 


For h = 0.1 and an iteration tolerance of € = 
io. the time per integration step varies between 15 
and 60 sec, due to the considerable number of itera- 
tions required. By setting « = 1078, the program is 
modified to correct only once, possibly at the expense 
of stability. The time per integration step is then 


reduced to 6 sec. 
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AUTONOMOUS DIFFERENTIAL EQUATION OF SECOND ORDER, 
FIRST DERIVATIVE ABSENT 


l. Purpose 


To integrate numerically a differential equation of 


the special form 
y" = fly) , (1) 
subject to the initial conditions 
= ‘ = 
y(%4) = Yo 7 Y'(%_) = 2y - 
Equations of the special form (1) occur, among other 
places, in problems of non-linear vibrations without 
damping. 
2. Method 
A simplified version of the Runge-Kutta method (see 


ENA, § 14.5). Let h be the integration step, and let 


x T= X + nh. Approximate values Yn of y (x) and ZH 


of y'tx) are computed by the formulas 
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hn? 
Yn+1 ~ Yn ba ena 2 fly) 
(2) 
Z =z + hf( + he ) 
n+l n Yn 2°n° ' 
n= 0, l, 2, ... . The error in a and am as h + 0 and 
x = x is fixed is o(h?). 


3. Flow Diagram 


We write x := Xi YrE Yo, 2S 2G y* denotes the 


current argument in f; A := hf(y* 
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h 
y + 52 


y + h(z + 3A) 


Display x 


Display y 
Stop 
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4. Storage and Program 


Ro R) Ry R 
h x y Zz 

00 

01 CLX 

02 STO 4 

03 RCL 2 

04 STO 7 

+ 05 \ 

06 

07 

08 

09 

10 Program to 

ll 

12 compute 

13 

14 £(y*) 

15 

16 

17 

18 

19 

20 RCL 0 

21 m 

22 STO 6 

23 RCL 4 


24 x # 0 


> 
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5. Operating Instructions 


Load the main program. Load the program for calcula- 
ting fly*) into program locations 05 + 19, assuming 
y* in Ro. Program must put F(y*) into X. Ro is availa- 
ble to store a constant needed for the calculation of 
f.. 

Move the operating switch to RUN. Select a mode of 
displaying numbers. If f contains trigonometric func- 


tions with arguments expressed in radians, press 


RAD 
Load data as follows: 
h into Ry 
Xo into Ry 
Yo into Ro 
29 into Ry 


Put auxiliary constants into R,- Pressing 
PRGM 
R/S 


will start calculation. The calculator will briefly 


display x1 and stop by displaying Yi- To exhibit zy 


press 


RCL 3 
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To continue computation, press 
R/S 


This will cause display of x, and Yo: etc. 


2 


6. Example and Timing 


We integrate the exact equation of the mathematical 


pendulum, 
y" = - 4 Ssiny , 


where g/2% = 9.81, under the initial conditions 


y(0) =F, y'(0) = 0 


(compare the program "Trapezoidal Rule"). The results 


for various values of h are as follows: 


Y 

x h=0.1 h=0.05 h=0.025 h=0.0125 
0.0 1.047198 1.047198 1.047198 1.047198 
0.1 1.004719 1.004785 1.004859 1.004885 
0.2 0.878363 0.879600 0.880033 0.880156 
0.3 0.675164 0.678906 0.680034 0.680340 
0.4 0.408987 0.416708 0.418874 0.419443 
0.5 0.102068 0.114850 0.118251 0.119124 
0.6 -0.216063 -0.198319 -0.193810 0.192679 
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Applying h?-extrapolation to the values at x 


yields 


0.102068 
0.114850 
0.118251 
0.119124 


0.119110 
0.119384 
0.119415 


= 0.5 


The theoretical hypotheses for continuing Romberg ex- 


trapolation beyond the second column of the Romberg 


scheme are not satisfied here. 


Computing time per integration step (including 


display of x) approximately 5 sec, 


choice of h. 


independently of 
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LINEAR SECOND ORDER DIFFERENTIAL EQUATION 


1. Purpose 


To integrate numerically the general linear second 


order equation, 


y" + p(x)y' + q(x)y = 0, (1) 
subject to initial conditions 
= t cs 1 
¥(Xq) Yo: y (xq) Yo: (2) 


2. Method 


Finite differences. Denoting by Y, an approximate va- 


lue of the exact solution y(x) at the point x, = Xp + 


nh, where h is the integration step, we approximate 
t 
y (x) by 


and 
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Introducing these approximations in the differential 


equation and writing 


Po. t= p(x) , Gd, te q(x) ; 


there results 


1 Ph = 
2 ¥na1 ay + yY Vet ThlYns1 ~¥ -1) + 4n¥, = Oe 


h n-l 


which can be solved for Yep t° yield 


= posi cf - 2 - _ih 
Yn+1 ~ (C2 A Cee ae (1 3P,)Yn-1! 
(3) 


1 


2 
= 2+ bp, “4 - 2h qn) Yn - (2 - hp) ¥,-3} . 


This is used to calculate the approximations Y, recur- 
sively. The starting value Yo = y (xq) is given. To ob- 


tain Y,, we use Taylor's expansion in the form 


2 


YR Yo. * Bg HO -¥6. 


Here Yo and Yo = ¥' (Xp) are given, and Yo can be cal- 


culated from the differential equation. There results 


2 
~ UL ce, he ‘ 
¥y = ¥o + hyg - 3-o¥9 + Yo] - (4) 


This calculation must be performed manually. 
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3. Flow Diagram 


Compute Yoel by (3) 


Display x, 


Display Yn 


Stop 
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4. Storage and Program 


Ro Bi Ro BS Ry Re Re Ry 

h x) Yo Yy temp temp 
xy Yn-1 Yn 
00 25 | 
01 } 26 «RCL O 
02 27 x? 
03 Program to 28 * 
a compute cd 
05 30 a 
06 p(x) 31 CHS 
07 32 4 
08 33 + 
09 RCL 0 34 RCL 3 
10 * 35 STO 
ll STO 4 36 7 
12 2 37 RCL 5 
13 RCL 4 38 - 
14 = 39 RCL 4 
15 RCL 2 40 2 
16 il 41 + 
17 STO 5 42 + 
18 43 RCL 
19 44 RCL 
20 Program to 45 + 
21 compute 46 STO l 
22 47 PAUSE 
23 a(x) 48 R4¢ 
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5. Operating Instructions 


Load the program, including the programs to compute 
p(x) (locations 01 through 08) and q(x) (locations 18 
through 25), assuming x in R,- The registers Re and R, 
are available to store constants needed in the calcu- 
lation of p and q. 

Move the operating switch to RUN. Select the mode 
of displaying numbers, and choose correct mode for ar- 
guments of trigonometric functions (ordinarily by 
pressing RAD) if such functions occur in p or q. Load 
data as follows: 


h into 
x) (:=x +h) into R 
Yq into Ry 


Calculate Yy from (4), and load it into R,- Pressing 


PRGM 
R/S 


starts computation. The calculator briefly displays Xo 


and stops while displaying Y>- Computation is conti- 
nued by pressing 


R/S 


which results in display of Xr Y3" etc. Continue in 
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this fashion until the desired range of x values has 
been covered. 

As always with differential equation problems, it 
h 


Nik 


is a good practice to repeat calculation with h := 
in order to check convergence of values Ys at fixed 
values of x. Romberg acceleration may be used at least 


once. 


6. Example and Timing 


We obtain values of the Bessel function Jy 00 by inte- 


grating the differential equation satisfied by it, 
l 1 
" = f) -— = 
y" + Ps y' + (1 >) y Oo, 
using 


_ = ES. apt _1l 
Yo = 540) =O, ys = T1(0) = 5 


No harm is done by the singularities of the coeffici- 
ent functions p and q at x = 0, because the program 
uses these functions only from x = x) onward. Equation 
(4), however, cannot be used to find Yy- From the po- 
wer series defining J, O09 we get to the same order of 


accuracy 
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The programs to calculate p and q are 


ol RCL 1 
02 1/x 
and 

18 1 
19 RCL 1 
20 xe 
21 1/x 
22 - 


fhe remaining locations are to be filled up with NOP 
instructions. (As an alternate possibility, the whole 
program could be compressed, which is not problematic 
here because there are no internal references. In this 
case, the instructions R/S and GTO 00 should be in- 
serted at the end.) 

Values obtained using several different steps h 


are as follows: 


Y J, (x) 

(exact 

x h=0.1 h=0.05 h=0.025 values) 
0.0 0.000000 0.000000 0.000000 0.000000 
0.1 0.050000 0.049958 0.049944 0.049938 
0.2 0.099667 0.099553 0.099517 0.099501 
0.3 0.148603 0.148406 0.148345 0.148319 
0.4 0.196436 0.196150 0.196063 0.196027 
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0.5 0.242806 0.242429 0.242315 0.242268 
0.287368 0.286899 0.286758 0.286701 

0.7 0.329790 0.329230 0.329063 0.328996 
3.6 0.095396 0.095459 0.095466 0.095466 
3.7 0.053631 0.053789 0.053824 0.053834 
0.012490 0.012740 0.012801 0.012821 

9 -0.027698 -0.027360 -0.027274 -0.027244 
4.0 -0.066613 ~0.066193 -0.066082 -0.066043 


Computing time per integration step approx. 4 sec. 


Part 6 


SPECIAL CONSTANTS 
AND FUNCTIONS 
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LOG-ARCSINE ALGORITHM 


l. Purpose 


To compute the functions Log x and Arcsin x for arbi- 
trary x > 0 and for arbitrary x ¢€ [-1, I], respective- 
ly, by an algorithm requiring only the evaluation of 


square roots. 


2. Method 


Numerical differentiation (see ACCA II, § 11.12). 
Clearly, for x > 0, 


Log x = o x" ” 
oh h=0 
Therefore, if 
wh) Ge ix song 
and 
8 eh (2 Py 40 nem OF Ly 2h. cas. y 
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then 

lim am = Log x . 

n-o 
Similarly, if x « (1, 1] and a := Arcsin x, then 

: oP ON an, eles 

Arcsin x = ah sin(ha) h=0 
Thus by putting 
sin(ha) 
a(h) := h 

and 

-n 

ay = a(2 ) , n=0, 1, 2, ... , 

we have 


lim an = Arcsin x . 


nro 


The sequences s_ = 2 and s_ = a_ both satisfy the 
n n n n 


same recurrence relation 


if 
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wo. Lye2 1 ak = <2 
Sai ia q (* =z) , So = 5 (x x) t (2) 
x 
then so = 2, and if 
Sy = xv¥l - an , S =x, (3) 
then s, =a. In both cases, the convergence of the 


sequence {s.} can be sped up by the Romberg algorithm. 
The foregoing algorithm fails for the trivial values 
Log 1 = 0 and Arcsin 0 = O, but it is very stable 


otherwise. 


3. Flow Diagram 


The program computes {20} if the index k is set equal 
to 0, and it computes {a)} if k # 0. In either case it 
saves the five most recent elements of the sequence 


{s}. 
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Pause and display So 
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4. Storage and Program 


0 1 2 3 4 5 6 7 
S_4 5.3 S_5 Si So x 
00 25 1/x 
Ol sTO 6 26 - 
02 RCL 7 27 4 
03 x= 28 + 
04 GTO 15 29 STO 3 
05 1 ~ 30 RCL1 
06 RCL 6 31 STO 0 
07 STO 4 32 RCL 2 
08 x? 33 sto l 
09 - 34 = =RCL 3 
10 vx 35 STO 2 
11 RCL 6 36 RCL 4 
12 * 37. «STO 3 
13. STO 3 38 2 
14 GTO 30 39 * 

+ 15 RCL 6 40 RCL 3 
16 ENTER 41 RCL 2 
17 1/x 42 + 
18 - 43 + 
19 2 44 vx 
20 = 45 RCL 3 
21 sto 4 46 * 
22 RCL 6 47 STO 4 
23 x? 48 PAUSE 
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5S. Operating Instructions 


Load the program, and switch to RUN. Select the mode 


of displaying numbers, such as 
SCI 8 


to get the floating eight-digit display. Load k into 
Ro, x into X (x #4 1 if k = 0 and x # 0 if k # OO). 


Pressing 


PRGM 
R/S 


will start computation. The calculator pauses briefly 
while displaying each sac No convergence test is pro- 
vided; convergence must be determined by inspection. 


Press 
R/S 


during pause to recover five most recent values. These 


will be in Ro through Ry, 


the Romberg algorithm. The most recent value is also 


the correct locations for 


displayed in the X register. 


6. Examples and Timing 


k := 0, x = 10 yields with FIX 9 display 


(2] 
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Log 10 = 2.302585094 in 16 iterations. 


Time required about 30 sec. There is an error of 
one unit in the last digit. By using SCI 8 dis- 
play we get 


Log 100 4.6051702 in 13 iterations, 
Log 1000 = 6.9077552 in 16 iterations. 


x 
i] 
oO 
“ 
i] 


e yields the elements S_yr Soe at oh 

$s, of the sequence {si} as follows: 

1.813430203 } 1.002606203 
1.175201194 Romberg 0.999991848 
1.042190612 acceleracion 1.000000050 

yields 

1.010449268 aad 1.000000001 
1.002606203 1.000000002 


Clearly, a FIX 8 display would have given the 


correct answer Log e = 1.00000000. 
k = 1, x := 1 yields with FIX 9 display 
Arcsin 1 = 1.570796319 in 15 iterations, 


which differs from the exact value 1/2 = 
1.570796327 by 8 * 107°. 
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THE GAMMA FUNCTION 


1. Purpose 


To compute the values of the gamma function, 


n“n! 
r(x) := lim = 
diag 
n-ew 
= fe’ +? at (xe > 0) 
0 


for arbitrary real values of x with a relative error 


Stirling's formula (see ACCA II, § 8.5). We select n 
such that Stirling's formula evaluates [(x+n) with the 


required low relative error and then obtain 


T (x+n) 


T(x) = 
(x) 


(see ACCA II, § 8.4). Stirling's formula is usually 


The Gamma Function 215 


given in the form 


es ee o(* - y¥2)Log x - x + J(x) : (1) 


where J(x) is the Binet function, 


Log ———- dt , 


to + x le- gers 


x > 0. It is known that J(x) possesses an enveloping 


asymptotic expansion as x > ®, 


ee eee ee 
x 360 x 1260 x 


I(x) % 


The direct implementation of Stirling's formula, 
approximating J(x) by the terms of the asymptotic se- 
ries shown above, would produce the required accuracy 
for x+tn > 5, but the resulting program would be too 
long for the HP-25. We therefore modify the formula, 


as follows. Writing (1) in the form 


ee ee /21. eX (Log x - £(x)) ; (3) 


we see that it is necessary to evaluate the function 


f(x) 1 - 2, (4) 
which within the prescribed tolerance may be approxi- 


mated by the polynomial in 1/x?, 
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g(x) :2 1 - — 5 | - —— (5) 
12 x 360 x 1260 x 


(see also Table 1). 

To evaluate g(x) directly by means of (5), using 
any of the available methods for computing the values 
of a polynomial, would require the storage of three 
constants with a total of 9 digits, occupying 9 pro- 
gram memory locations. The resulting program would 
again exceed the capacity of the HP-25. 

It turns out that a shorter program is obtained by 
converting the function g(x) into a continued fraction 
(i.e., by representing it by its diagonal Padé appro- 
ximant), using the qd algorithm (see ACCA, § 12.4). 
The following qd table results from the four coeffi- 
cients of the function g: 


(n) (n) (n) 
a ay ey G2 
at 
oe 12 - 
12 a 20 ea 
1 30 _ 53 315 
360 2 210 
~ 7 
1260 


This yields the continued fraction 
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x 
deca ‘ 
ee 
= 
os f20 
23 
315 
5 
x 


which agrees with g(x) up to an error of at most 


o(x 8) as x + », In the above fraction we may replace 
53 53 


the unwieldy 315 by 318 7 - without affecting the ac- 


curacy of the Final result. ( would arise exactly if 
a in g(x) were replaced by eae ~ 
l 1260 21600 

T2790") By algebraic manipulation the resulting appro- 


the coefficient 


ximation to £(x), 


x 
h(x) := oe ioe Ne 

x + ke = 

_ _20 
* i 
6 

yaaa 
x 
May be simplified as follows: 
h(x) = + ________-— : 
1+ T 
3x(4x + ) 
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This, finally, is the approximation to f(x) used in 
our program. Only 4 one-digit constants need to be 
stored. The numerical values given in Table 1 show 
that the error committed in approximating g(x) by 
h(x) for x 2 5 again lies in the permitted range. If 
x < 5, our program thus will determine the integer n 


< 
mentioned initially such that 5 = xtn < 6. 


x £(x) g(x) h(x) 

1 0.918938533 0.918650794 0.918566776 
2 0.979329652 0.979327877 0.979327701 
3 0.990774025 0.990773946 0.990773914 
4 0.994802332 0.994802324 0.994802324 
5 0.996671061 0.996671060 0.996671061 
6 0.997687312 0.997687312 0.997687312 
7 0 .998300470 0.998300470 0.998300470 
8 0.998698592 0.998698592 0.998698592 

Table 1 


3. Flow Diagram 


We let y := xtn, p := (x). 
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ry) i= /# 


« ey (Log y ~ hly)) 
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4. Storage and Program 


R R R R R 


0 1 2 3 4 
x P Y 
00 
ol STO 0 
02 1 
03 STO 1 
04 5 
05 RCL 0 
> 


08 STO*1 
09 1 
10 + 
11 GTO 06 
> 12 STO 2 
13 6 
14 * 
15 1/x 
16 RCL 2 
17 7 
18 5 
19 * 
20 1l/x 
21 RCL 2 
22 4 
23 * 
24 + 
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5. Operating Instructions 


Load the program; operating switch to RUN. Select a 


mode of displaying numbers, ordinarily by pressing 


SCI 8 


to get 8 digit floating display. Load x into X regi- 


ster. Pressing 


PRGM 
R/S 


will cause the value of [(x) to be displayed. Values 
are generally correct to 8 significant digits; for 
large x errors may be slightly larger due to inaccu- 
racies in exponential function. 

To get [(x) for new value of x, load new x into X. 


Pressing 


R/S 


will yield new value. 
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6. 


(4) 


Special Functions 


Examples and Timing 


yields r(x) = 1.0000000 
1.0000000 
2.0000000 
6.0000000 
2.4000000 * 10 
1.2000000 * 10 
7.1999999 * 10 
5.0400000 * 10 
4.0319999 * 10 
3 


-6287999 * 10 


wo mon Dn WW BP WwW NY 


nN & WN NY FF 


= 
Oo 


These values provide a genuine check, because the 
relation [(n+l) = n! is not used in the program. 


Computing time for each value about 5 sec. 


x = 0.5 yields rg) = 1.7724538, which agrees 
with the exact value, /7. 


x = 0 correctly yields "Error", because [(x) has 
a pole at x = 0. The same holds for all negative 


integers. 


x = -2.37 yields [(-2.37) = 1.2183595. 
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INCOMPLETE GAMMA FUNCTION 


1. Purpose 


To compute, for arbitrary x > 0 and arbitrary real a, 


the values of the incomplete gamma function, 


P(a,x) r= st? “eo dt. (1) 
x 


2. Method 


Expansion in terms of reciprocal Laguerre polynomials. 


Let 


denote the Laguerre polynomial of order n with para- 
meter a (see ACCA I, p. 119). Then there holds for all 


x > 0 and alla<il 


[(a,x) = x7 e * t snus 1 (2) 
i 7 “ (ntl)! . (-a) )_ (-a) ,_ 
n=0 Ln (-x) Lied (-x) 
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(see ACCA II, § 12.13). This expansion also holds for 


a 2 l, provided -x is not a zero of one of the polyno- 
(-a) , 
n ’ 
exceptional case. Using the asymptotic theory of the 


mials L see example [4] for how to deal with the 
Laguerre polynomials (see Szegd, Orthogonal polynomi- 
als, Theorem 8.22.5) the n-th term of the above series 
(including the factor xe %) can be shown to be asym- 


ptotic to 
an /® 7 4vnx 
a ‘ 


The convergence of the series (2) thus is subgeometric. 
More precisely, the number n of terms required to get 


k correct decimal digits is asymptotically equal to 


n= 0.33 — (3) 


and thus is roughly proportional to K?, independent of 
a, and inversely proportional to x. 
For purposes of computation we write (2) in the 


form 


I(a,x) = £ a ; (4) 
n=1 “n-l'n 
where 
eae 
Yn is x76 755 n-l 
n-1l 


and 
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a_-x _. Roa (5) 


From the recurrence relation of the Laguerre polynomi- 
als (see Szegé, Orthogonal polynomials, equation 
(5.1.10)) we find 


}, 


i ars (n-a- ljhs 


9. a (Sy Se st ee 
n n 
n= 1, 2, ... . To save arithmetical operations, this 


is used in the form 


ae _ 
£8 = al (n -a- 1) (2, & _») + (n+ x)£ a} (6) 


-1 
Summation is stopped if 


et 
anol 


where e€ may be chosen arbitrarily. 
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3. Flow Diagram 


Display 
T'(a,x) 


To save operations, most operations indicated in the 


n := n + 1 box are anticipated in the box marked *. 
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4. Storage and Program 


0 1 2 3 4 5 6 7 
42 aad ss ‘nh : [e] 

00 25 RCL 6 
Ol CLX 26 RCL 1 
02 sto 2 27 + 
03 sto 4 28 RCL 
04 1 29 sto 2 
05 Sto 3 30 * 
06 STO 6 31 ‘f 
07 RCL 1 32 RCL 6 
08 RCL O 33 : 
09 y~ 34. STO 3 
10 sto 5 35 RCL 2 
Li’ Reis f 36 * 
12 e* 37 + 
13. STO5 38 «= STO+4 

>» 14 RCL 5 39 PAUSE 
15 RCL 6 40 ABS 
16 RCL 0 41 RCL 7 
i? - 42 xey 
18 sTO*s 43 GTO 49 
19 1 44 1 
20 ~ 45  STOo+6 
21 RCL 3 46 RCL 6 
22 RCL 2 47 STOs5 
23 - 48 GTO 14 


24 i + 49 RCL 4 
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5. Operating Instructions 


Load the program, move the operating switch to RUN. 
Select a mode of displaying results, for instance by 


pressing 


FIX 
9 
Load data as follows: 
a into Ro 
x into Ry 
€ into Ry 
ae -10 F 
For 9 digit accuracy, let e« = 10 . Pressing 
PRGM 
R/S 


will start computation. The calculator briefly dis- 
plays each term of the series (4) and stops if that 
term is < ¢€, displaying value of [(a,x). To get new 


value, it suffices to load new data and to press 


R/S 
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6. Examples and Timing 


The incomplete gamma function is of interest because 
it contains as special cases a number of important 

functions of mathematical physics and of statistics. 
The following special cases occur in the Handbook of 


Mathematical Functions by Abramowitz and Stegun. 


The exponential integral, 


obviously is given by 


E, (x) = T(0,x) . 


Using the tolerance ec = rete the Following va- 


lues are obtained for different values of x: 


number n computing 
of terms time 


x E, (x) required in sec 
5 0.559773594 67 190 
: 0.219383934 36 100 
0.048900511 19 54 
4.0 0 .003779352 10 27 


The computed values are in error by at most one 


unit in the last place. 
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The exponential integrals of higher order, 


wo e *t 
E(x) := f nr at, neal, 2, % 
¥ lit 
may be expressed as 
E(x) = xl p(i-n,x). 


Some sample values computed by our program (using 


€ = 107?) are as follows: 


x E, (x) E,¢%) E59 () 
0.3266439 0.1652428 0.0310612 
1.0 0.1484955 0.0860625 0.0183460 
0 0.0375343 0.0250228 0.0064143 
0.0009965 0.0007830 0.0002783 


All values are correct to the number of digits 


given. Computing times are comparable to those 
given under 2 


The error function, 


> x 
erf(x) = —fe dt, 
yn 0 


equals 
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as may be seen by a simple substitution. Using 


the tolerance « = aaa the program permits the 


following values to be obtained: 


number n computing 
of terms time 


x erf(x) used in sec 
1. 0.842700793 34 97 
2ni 0.995322265 10 30 
3.0 0.999977910 5 15 


The values given are correct; however, the pro- 
gram Error function is much more efficient for 


x << 2, 


Chi-square probability function. This function 


may be defined as 


2 
; e 
Q(x" ]v) r= 
) 


Vv 
(3 


Some values provided by our procram (using € = 


10°") are as follows: 


v x? 

1 2 3 4 
2 0.60653 0.36788 0.22313 0.13534 
4 0.90980 0.73578* 0.55783 0.40601 
6 0.98561 0.91970 0.80885 0.67669* 
8 0.99825 0.98101 0.93436 0.85712 
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As mentioned in section 2, for each a 2 l the 
program may fail for isolated values of x due to 
negative zeros of the Laguerre polynomials. For 
instance, in the above table (x2 ] v) = (2|4) is 
such a pair of values. Failure to compute [ (a,x) 
at x = Xo is recognized by an error halt. In such 
cases, P(a,Xp) may be obtained as the arithmetic 
mean of the values Pla, xX9*5) where 6 is suffi- 
ciently small. In the above table, the values 
marked * were obtained in this manner. Of course, 
near failure may occur for values near XQ: This 
is recognized by the appearance of large terms in 
the series (4). Again, more precise values can 


then be obtained by linear interpolation. 
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ERROR FUNCTION 


L; Purpose 


To compute the error function, 


x 
fe at, (1) 
0 


-9 
for all x 2 0 with an error of less than 10 


2. Method 


A power series is used for 0 x < 2, and a continue 


> 
fraction for x = 2. 


a) The Taylor series x = 0, 


ro) 
z 


es (-1)"_ antl on 
Ya n=0 


aoe ni (anely 


is not to be recommended for numerical computation, be- 
cause the large alternating terms may cause cancella- 
tion, and because the general term is awkward to com- 


pute recursively. However, the series (2) may be ex- 
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pressed as a confluent hypergeometric series (ACCA I, 
§ 1.5), 


erf (x) 


a 
a| 


We thus may use Kummer's first identity (ACCA I, eq. 
(1.5-3); ACCA II, eq. (9.7-8)) to obtain 


= 2% . 3, 
erf(x) = e 1F1 (l ; 5 x ye ok (3) 


In the last series, all terms are positive, and no 


cancellation can occur. Moreover, if 


2 7 
= e* 1Fy (l ; 3 ; “7) = aye 
YT n=0 
then 
oe anigen 
0 
2 
a= x a 
n 1 “n-1 ’ 
n+ 5 


n=1, 2, ... . This formula is used in the program to 
evaluate the series (3) for 0 =x < 2. 
b) Numerical experiments show that the series (3) 


can be used safely to compute erf(x) at least up to 
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x = 5.0, where erf(x) = 1 to the accuracy required. 
However, the time required to compute erf(x) by (3) is 
about 10x seconds for x 2 2 due to the increasing num- 
ber of terms that have to be summed. For x 22 we 
therefore use an alternate method that is considerably 


faster. For all x > 0 there holds 
erf(x) = 1 - erfc(x) , 


where 


is the complementary error function. This function ad- 


mits the asymptotic expansion 


which however cannot be used to obtain arbitrarily ac- 
curate values for any x. An expression that converges 
for all x > 0 is obtained by converting the asymptotic 
series into a continued fraction (ACCA II, § 12.13). 


This yields 


‘ > 
It has been determined experimentally that for x = 2 
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the 16th approximant yields erfc(x) with an error < 
-9 


10 ~. Thus our program uses the approximation 
oe 1 | 
erf(x) = 1-£ aa , + 2 + + 8.) 

ae Gee fe 


Working with a fixed number of approximants has the 
advantage that the continued fraction can be evaluated 


rapidly working from the tail upward (ACCA II, § 12.1). 


3. Flow Diagram 


e = Sef 1 zg - i 
We set a := ave s:=f ane vr=nt 5, mrs 8 5k. 


Error Function 


Display 


erf (x) 


Display 


erf (x) 
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4. Storage and Program 


Ry R R R R R R R 


1 2 5 4 5 6 7 
[x| a v s m aa 5 E 

00 25. ~STO*l 
01 STO 0 26 RCL 1 
02 x? 27. STO+3 
03 CHS 28 RCL 7 
04 e* 29 x<y 
05 RCL 5 30 GTO 18 
06 : 31 RCL 3 
07 «STO l 32. GTO 00 
08 2 » 33 RCL O 
09 RCL O 34 8 
10 xZy » 35 sto 4 
1l GTO 33 36 x Sy 
12 sTO*L 37 : 
13. RCL 6 38 RCL O 
14 s$To 2 39 + 
15 STO#1 40 RCL 4 
16 RCL 1 41 RCL 6 
17 STO 3 42 a 

+18 RCL O 43 x #0 
19 x? 44 GTO 35 
20 RCL 2 45 R + 
21 1 46 sTorl 
22 + 47 1 
23. STO 2 48 RCL 1 
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5. Operating Instructions 


Load the program; move the operating switch to RUN. 


Press 


FIX 


for nine-digit display (results will be accurate to 9 


digits). Load € := 162° into Ro: 


Load 0.5 into Re: 


STO 6 


Load Ym into Re? 


sTo 5 
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> 


To obtain erf(x) for arbitrary x 0, load x into 


X register and press 


PRGM 
R/S 


Value of erf(x) will be displayed after a few seconds. 


To obtain erf(x) for new x, load new x into X and press 


R/S 


6. Examples and Timing 


= 1.99 yields erf(x) 
= 2.00 yields erf(x) 
0.01 yields erf (x) 
= 5.00 yields erf(x) 


0.995111413, time 20sec. 
0.995322265, time 11sec. 
0.011283416, time 3 sec. 
1.000000000, time 11sec. 


The changeover point from power series to conti- 


x Re KR 
iT} 
“ou 


Mea 


nued fraction may be changed to any integer k 
such that 0 < k < 10, by substituting k for 2 in 
instruction 08. It is thus possible to compare 
the performance of the formulas (3) and (5) for 
various values of x. For instance, for x = 4 we 


get 


by (3) erf (x) 0.999999985 , time 40 sec, 
by (5) erf(x) = 0.999999985 , time 1] sec. 
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COMPLETE ELLIPTIC INTEGRALS 


1. Purpose 


To compute, for 0 = k < 1, the complete elliptic inte- 


gral of the first kind, 


n/2 
K(k) := ff —————_ do, 


2. Method 


The arithmetic-geometric mean. Let a), By be positive 


numbers, and let 


n= 0, l, 2, ... . Then the sequences {a} anda {80} 
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rapidly converge to a common limit u(ag,BQ), called 


the arithmetic-geometric mean of XG and Bo: If A = 1, 
Bo =k' := Yl - Ke; then 
T 

K(k) = Mk) 

Moreover, if Yo °= k, 
meee 
Vuel or 5 fo, 8.) , 
n= 0, l, 2, , then 
B(k) = K(k)(L- 5 2™y) 
m=0 


(see ACCA II, § 9.10, problems 9 - 13). 


3. Flow Diagram 


We let a :=a_,8:= 8, a' 


n 
k_2 
I 2 Yr T:= 2 
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Yes 


244 


4. 


Special Functions 


Storage and Program 


os) 

a 

Ct 
wn IN 


RCL 
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5. Operating Instructions 


Load the program; move the operating switch to RUN. 
Select a mode of displaying numbers, for instance by 


pressing 
FIX 9 


2 
to get 9 digits after decimal point. Load «€ = 10 


into R as follows: 


1’ 


STO 1 


Load k into Ro and also leave it in X register. Pres- 


sing 


PRGM 
R/S 


will after a few seconds cause K(k) to be displayed. 


To display E(k), press 


AV 


K(k) then will be in Y register. 
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6. Examples and Timing 
For k = 1/72 = 0.707106781 we get 


K = 1.854074677 , E = 1.350643881 ; 


computing time approx. 8 sec. For k = 0.9999 the re- 
sults are 


K = 5.645148167 , E = 1.000514502 


; 


computing time approx. 12 sec. 
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BESSEL FUNCTIONS, INTEGER ORDER 


Ls Purpose 


To compute the Bessel functions J tx) for arbitrary 


x > 0 and for the orders n = 0, 1, 2, 3, 4. 


2. Method 


J.C.P. Miller's method of backward recurrence (see W. 
Gautschi, Computational aspects of three-term recur- 
rence relations, SIAM Review 9, 24 - 82 (1967)]. We 
use a recurrence relation satisfied by the sequence 


to Ce) } (x fixed) in the form 


2n 
= “fn = (1) 
Jn-1 x oF Jnel 


If the recurrence 
2n (2) 


Yn-1 ~ &% Yn ~ Ynt1 


is solved with arbitrary starting values, Say, 
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where the integer m is sufficiently large, then the 
elements of the sequence y 


J 


».+.- SOON become 
m-n’ Ym-2? 


proportional to J -.. , because {Jo} is a 


m-1' “m-2" 
"recessive" solution of the difference equation (2). 


Thus, approximately, 
Yor (3) 


where c is independent of n. The factor c may be de- 


termined from the known relation 
JQ (x) + 23, (x) + 25, (x) + ...21, 
which gives 
C= Y¥_ + 2y, + 2y, + --- - (4) 


In this program m may be selected as an arbitrary 
positive integer; see the examples for a rule of thumb 
on how to select m as a function of x. The program 
then constructs the sequence ty}, always saving the 
last five elements. The values JOO) are then computed 


from (3), approximating c by 


+ cee . 
Yo 2y> + 2¥4 + + 2¥2 (n/ 2] 
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3. Flow Diagram 
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4. Storage and Program 


Ro Ry am Ry Ra Be Re 
oF J, os a5 Jy c x 
00 25 2 
01 STO 6 26 . 
02 2 27 «sTo 7 
03 sTo*7 28 x = 0 
04 0 29. GTO 40 
05 sTol 30 4 
06 STO 5 31 
07 1 32 FRAC 
08 sto 0 33 x #0 
+09 RCL 3 34 GTO 09 
10 STO 4 35. RCL O 
ll RCL 2 36 2 
12 $T0 3 37 * 
13. RCL 1 38 sTOrs 
14 sT0 2 39 GTO 09 
15 CHS >» 40 RCL O 
16 RCL O 41 STO+5 
17 sTo l 42. RCL 5 
18 RCL? 43 1/x 
19 RCL 6 44  sTo*0 
20 ¢ 45. sTO*l 
21 * 46 sTO*2 
22 + 47 §TO*3 
23 sto 0 48 stor 


2n 
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5. Operating Instructions 


Load the program; turn the operating switch to RUN. 
Select the mode of displaying numbers, for example, 


by pressing 


SCI 8 
Load m into Ro. (m = 30 is adequate for values x < 10.) 
Load x into X. Press 

PRGM 

R/S 


to start the computation. The calculator will stop by 


displaying calculated value of J9(x); the values Jy, (x) 
(k = 0, 1, ... , 4) are also found in R,. To restart 
the computation with a new value of x, simply reload 
mM into Ro, the new x into X, and press 


R/S 


6. Examples and Timing 
For x = 1, m = 30 we get the values 


0.76519769 
0.44005059 
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Jy = 0.11490348 
J, = 0.01956335 
J, = 0.00247664 


which are exact to the number of digits given. 


Computing time 40 sec. 
For x = 10, m= 30 the results are 


= -0.24593576 
= 0.04347275 
= 0.25463031 
= 0.05837938 
= -0.21960269 


correct to all digits given. 


For x = 17.5 we choose m = 40. The computing time 
now is about 59 sec. There results 


= -0.10311040 
= -0.16341997 
= 0.08443383 
= 0.18271913 
= -0.02178727 


QQ GO 4G 
hb WN fF Oo 


Again, the errors are < Lo The value m = 50, 


which theoretically produces still more accurate 


results, causes overflow. 
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BESSEL FUNCTIONS, ARBITRARY ORDER 


1. Purpose 


To compute, with an error < 10 8, the values of the 


Bessel function J (x) of arbitrary real order v for 
< < 
0 =x = 10. 


2. Method 


Power series expansion. We have (see ACCA II, § 9.7) 


(e1yPeo? 
4 


(v+l) ont ; 


x,V 
7 (5) 
~ ~P(vtl) 


a 
z 


J. (x) 
“ n=0 


We write this in the forma f= Boe where 


Gy” 
= Tory 2 a= 


and the remaining bo are calculated recursively by 


x, 2 
(5) 


A open) ned 
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The summation is terminated if |b | < 1071? phe term 


r'(v+l) required in the computation of a must be compu- 
ted separately, if necessary by the program gamma 
function. 
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3. Flow Diagram 


We write s := 2 b 


| Display J, (x) = 2» 


256 Special Functions 


4. Storage and Program 


R R R R R R R 


0 1 2 3 4 5 6 

ig OG 
00 25 RCL 6 
01 2 26 * 
02 < 27 + 
03 STO 7 28 CHS 
04 1 29 =sTo*S 
05 EEX 30 RCL S 
06 CHS 31 SsTOo+4 
07 1 32 ABS 
08 2 33. RCL 2 
09 STO 2 34. x2 y 
10 RCL 7 35. GTO 39 
11 RCLO 36 1 
12 y* 37 STO+6 
13. RCL 1 38 GTO 20 
14 + + 39 RCL 3 
15 sTo 3 40 RCL 4 
16 1 41 * 
17 «sT0 4 42 GTO 00 
18 sTO 5 43 
19 STO 6 44 

+20 RCL 7 45 

21 x? 46 
22. RCL 6 47 
23 RCL O 48 


24 + 49 


NIX 
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5. Operating Instructions 


Load the program, and move the switch to RUN. Press 
FIX 8 
to get a reliable eight-digit output. Load 


v into Ro 


T(v+l) into Ry je 


(The latter value may be computed by the program gamma 


function.) Load x into the X register. When 


PRGM 
R/S 


is pressed, the calculation will start. The calculator 


stops by displaying J (x) : 
This program does not work if v is a negative in~ 


teger, v = -n. Here the relation 
(x) = (-1) "7. (x) 
TS n 


should be used. The program may be inaccurate if v is 


close to a negative integer. 
For nonintegral v, the Bessel function of the se- 


cond kind may be computed using the relation 
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J (x) cosvt - J (x) 
v -v 
Yo (x) = —— i 
v sinvt 
For integer values of v,v =n, the limit of the fore- 


going as v + n may be calculated numerically, using 
the Romberg algorithm (see example [4]). 
6. Examples and Timing 


Jo (2) = 0.22389078, all digits correct, time re- 
quired about 13 sec. 


J, (10) = -0.24593576, all digits correct, time 
required about 27 sec. 


To compute J,,5(5), we require 


5 1 3 
ry 5 +5 0 
Computation yields F426) = ~0.16965131, in 


agreement with the exact value. 
[4] To compute ¥) (x), we may use the relation 


-2 2 
Yq (x) = av J Ok) v=0 ° 


Letting 
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Our program yields for x = 2 


Vv qd, 

2 0 .00000000 
1 0.57672481 
0.5 0.74780184 
0.25 0.78844828 
0.125 0.79839964 


With the Romberg algorithm the limit of these va- 
lues as v >» 0 is determined as 0.80169622. Thus 


we obtain 
Yo (2) = 0.51037566 , 


which is in error by only l * Mig 
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BESSEL FUNCTIONS: ASYMPTOTIC SERIES 


1. Purpose 


To compute the asymptotic series, for arbitrary x > 0, 
of the Bessel functions of both the first and second 
kind, 


J (x) and YX (x) ' 


for arbitrary real orders v. 


2. .Method 


The asymptotic series is computed most conveniently 
from the formula 


. qT 
i(x - =9) 
F u /2 2 #4 1 lL 1 
Tyo) +i ¥o0 */e Boe as eggs te 


where is a hypergeometric series (see ACCA II, 


2Fo 
§ 11.5). We write 
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and calculate the ay by means of the recurrence 


(- ; +v +n)(- 5 - v +n) 


o¢ n? 2xn 


The sums R and § are then given by 


The point of termination of these series may be chosen 
arbitrarily. Ordinarily, to get best results, one 


should terminate as soon as lal > |a If v equals 


nail 
an integer plus y¥2, the sums R and S will terminate 
mathematically, and the representation given by (1) is 
not only asymptotic but exact. 

To save programming storage, the correct assign- 
ment of the terms ta, to the sums R and S is left to 
the operator (see operating instructions). It then be- 
comes possible to compress the whole program into 49 


steps by converting R + iS into polar coordinates, 
R+is=te'*® (r>0), 
which yields 


i(x - 
e (R + iS) =Te 
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Converting back to cartesian coordinates then yields 


Nu 2 vial T 

J (x) % ee T cos(x - > a o) , (2) 
% 2 ; va T 

YO) ae od T sin(x - ener + >) . (3) 


3. Flow Diagram 


We set q := (-Sevtn) (-g-ven) /2xn. 
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n+l 


Display q, a := ga 


Summation to 


be continued? 


manually 


n even? 


(n+1)/2, 


S§ :=§ + (-1) 


Display (2), (3) 


264 


4. 


> 


Special Functions 


Storage and Program 


PAUSE 


> 
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5. Operating Instructions 


Load the program, and switch to RUN. Press 
FIX 8 


to obtain adequate display of results. Press 


to cause arguments of trigonometric functions to be 
interpreted in radians. Load 


x into Ry 


2v+l into Ry 


(the index v must be stored in this form for reasons 


of space). Pressing 


PRGM 
R/S 


will first cause qn to be computed and to be displayed 

briefly. If |q_| < 1, allow computation to continue by 
n 

doing nothing. The machine will stop by displaying an: 


Press 


4k+1 
4k+2 


sto-5 if n 
sTo-4 ifn 
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STO+5 if n = 4k+3 
4k+4 


STO+4 ifn 
Pressing 
R/S 


will cause the computer to proceed and to calculate 


next q_ and next a. If la, 2 1, quickly press 


R/S 
GTO 29 
R/S 


At next stop, the calculator will display asymptotic 
value of J (x). To get asymptotic value of YOO, 


press 


AV 


Caution: If R < 0 and |S/R| < Lor oy the value of 9 


will be in error by precisely 1, and signs of both Jy 


and xy will be wrong. 


6. Examples and Timing 


v = 0, x = 10. Stopping at n = 20 yields 


J, (10) % -0.24593576 , ¥Q (10) ®% 0.05567117 
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Both values agree with the exact values to all 


digits. Total computing time about 95 sec. 


J,5,(2) = 0.51301614 , Yy2 2) = 0.23478571 


These are the exact values. 


- il = 
MBS kD 


Jyy72 69) = 0.19056437 , Yaas2'>) = -0.57174942, 


These values agree with those computed by the 


program Bessel functions, arbitrary order. 
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RIEMANN ZETA FUNCTION ON CRITICAL LINE 


1. Purpose 


To compute the function 


: erGe a a 
S(t) r= n(5 + it) 


for real values of t, where 


nis) r= SEU) (8) n 7 gis) , 


and ¢ denotes the Riemann zeta function, 


G(s) := I on (Re s > 1) 


The function € is entire, even, and real for t real. 


There holds €(t) = 0 if and only if s = 5 + it isa 


nontrivial zero of ¢(s) (see ACCA II, § 10.8). 


2. Method 


We use the integral representation 
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&(t) = Sf £(x) cos(2tx) ax , (1) 
0 
where £(x) := $"(x) - (x), 
g(x) r= e% wlel®) , 
) 2 
w(u) := ZI oe Ta 
n=l 


(see Bieberbach, Funktionentheorie Vol. II, Chapter 8). 
The integral (1) is approximated by the midpoint 


formula, using an arbitrarily chosen step h, that is, 


(2) 
E(t) =h £ f(x, ) cos (2tx, ) ; (2) 
k=0 


where x TF 


sh + kh. The summation is terminated as 
soon as | £(x,) | < €. 


The function f is evaluated by 
means of the series 


% e » (3) 
1 


The summation is stopped when the general term (which 
is always positive) becomes < €. The series converges 
rapidly for all x = 0. 

Because of lack of programming space, no automatic 


step refinement is built into the program. However, 
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the step should always be refined manually to check 
the accuracy of values. Because f(x) = O(exp ce 4%) 
as x > ©, the error of the numerical values of the 
integral (1) is o(h™) for every m > 0, and Romberg 


acceleration is therefore neither required nor effec- 
tive. 


3. Flow Diagram 


We denote the sum in (2) by £, the sum in (3) by S, 


and the n-th term of S by au: 
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n+l 


compute an 


Xo:= f + £(x)cos(2tx) 


Display €(t) 
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4. Storage and Program 


RCL 3 


- 
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5. Operating Instructions 


Load the program, and switch to RUN. Press 


SCI 8 


to get the floating eight-digit display of results. 
Load the data as follows: 


2t into 


Ro 
h into Ry 
ce into R, 
5h into R3 . 
Clear R, by 
CLX 
STO 4 


Because argument of cosine is in radians, press 
RAD 


When 


PRGM 
R/S 


is pressed, the calculator will start whirling and 
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stop by displaying €(t) as computed by formula (2). To 


check accuracy, the computation should be repeated 
1 


with h := 5h. The necessary new data are furnished by 
RCL 1 2 
2 + 
+ STO 3 
STO 1 CLX 
STO 4 


Note: Some of the foregoing instructions can be 
automated on a calculator with a few more memory lo- 


cations. 


6. Examples and Timing 


Using € := 10° , we get for 


t = 14 h=0.2 — € = -0.46487914 ( 24 sec) 
0.1 +0 .00063699 ( 49 sec) 
0.05 0.00020129 ( 99 sec) 
0.025 0.00020129 (198 sec) 
t = 15 h=0.1 — € = -0.00001633 
0.05 -0.00070570 
0.025 -0.00070570 


Thus a zero exists between t = 14 and t = 15, implying 
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the existence of a zero of the Riemann zeta function 
on the critical line between s = 5 + 14i and s = ; + 
15i. 


INDEX 


Aitken acceleration, 48, 53 
asymptotic expansion, 162, 235 


for Bessel function, 260 


Bell number, 150 
Bernoulli number, 132 
Bernoulli's method, 80, 87 
Bessel function, 
of arbitrary order, 253, 260 
of integer order, 201, 247 
of second kind, 257, 260 
binary expansion, 126 


binomial coefficient, 18, 72 


chi-square probability function, 231 
continued fraction, 24, 31, 235 


convergence, quadratic, 53 


deflation, 74, 77 
a@ifference, finite, 196 
differential equation, 
autonomous, of second order, 189 
of first order, 181 
linear, of second order, 189 
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elliptic integral, complete, 241 
Enestr6ém-Kakeya theorem, 118 

equation, 41 

error function, 230, 233 

Euclidean algorithm, 14, 19, 33, 34, 36 
Euler's constant, 167, 175 

exponential integral, 229 


of order n, 230 


fixed point, 41, 53 


Fourier transform, fast, 124 


gamma function, 214 
incomplete, 223 
Gauss-Lucas theorem, 113 


greatest common divisor, 17, 19 


Horner algorithm, 67, 73 
group property of, 72 


integration, numerical, 155, 166 
irrationality, quadratic, 30, 31 
iteration, 41, 47, 53, 182 


Lagrange-Btirmann theorem, 140 
Laguerre polynomial, 223 


Log-Arcsine algorithm, 207 


map, contracting, 42 
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mean, arithmetic-geometric, 241 
midpoint formula, 155, 269 
Miller, J.C.P., formula of, 135 


Newton's method, 
for complex roots, 59 
for polynomials, 73 


for reciprocal power series, 124 


pendulum, 186, 194 
Plana summation formula, 170 
polynomial, stable, 101 

with zeros in unit disk, 107, 114 
power series, 

exponentiation of, 144 

power of, 134 

reciprocal of, 123 

reversion of, 140 


prime factor, 9 


quotient-difference algorithm, 
87, 95, 216 


recurrence, backward, 247 

Riemann zeta function, 177 
on critical line, 268 
zeros of, 275 

Romberg algorithm, 156, 160, 161, 162, 169, 
171, 176, 177, 178, 186, 188, 195, 209, 259 
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Routh algorithm, 101 
Runge-Kutta method, 182, 189 


Schur-Cohn algorithm, 107, 114 
Stirling's formula, 142, 214 


trapezoidal rule, 181 


The book contains some 35 programs written for a specific programmable 
pocket calculator, the HP-25. These programs implement algorithms in 
number theory, equation solving, algebraic stability theory, calculus of 
power series, numerical integration as well as algorithms for the evalua- 
tion of special higher transcendental functions, such as the Gamma func- 
tion, various Bessel functions, the error integral, and the Riemann zeta 
function. 

A unified format has been chosen for the preparation of the programs. 
By means of the flow diagrams and the detailed descriptions that are 
provided, the programs are easily adapted to run on any calculator of 
comparable capacity. 

The main benefit for users of these programs consists in the immediacy 
and concreteness in which they will experience the performance of the 
numerical algorithms presented. Convergence, speed of convergence, 
accuracy, and reliability of an algorithm will be experienced much more 
vividly when the result becomes visible at once on the users’ own calcu- 
lators. This is the first collection of high-level mathematical programs for 
a pocket calculator. 
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